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ABSTRACT 

As  a  first  step  to  understand  the  compressibility  effects,  interaction  of  isotropic 
quasi-incompressible  turbulence  with  a  weak  shock  wave  was  studied  by  three- 
dimensional  time-dependent  direct  numerical  simulations.  In  addition,  linear  anal¬ 
ysis  was  used  to  study  interaction  of  isotropic  turbulence  with  shock  waves  of  a 
wide  range  of  strengths.  The  effects  of  the  fluctuation  Mach  number  A/<  and  the 
average  Mach  number  of  the  upstream  turbulence  on  turbulence  statistics  were 
investigated. 

Both  numerical  simulations  and  linear  analyses  of  the  interaction  show  that  tur¬ 
bulence  is  enhanced  during  the  interaction  with  a  shock  wave.  Turbulent  kinetic 
energy  (TKE)  and  transverse  vorticity  components  are  amplified,  and  turbulent 
length  scales  are  decreased.  The  predictions  of  the  linear  analyses  compare  favor¬ 
ably  with  simulation  results  for  flows  with  Mt  <  —  1,  which  suggests  that  the 

amplification  mechanism  is  mainly  linear. 

Rapid  evolution  of  TKE  just  downstream  of  the  shock  was  not,  however,  repro¬ 
duced  by  the  linear  analysi<^  Investigation  of  the  budget  of  the  TKE  transport 
equation  shows  that  this  benavior  of  TKE  is  manifested  in  the  pressure  transport 
term  which  is  nonlinear.  The  budgets  of  enstrophy  components  show 

that  their  amplifications  through  the  shock  are  mainly  caused  by  the  distortion  due 
to  the  mean  flow  compression,  and  that  effect  of  baroclinic  torque  is  not  significant. 

Shock  waves  were  found  to  be  distorted  by  the  upstream  turbulence,  but  still 
have  a  well-defined  shock  front  for  Mi  <  - 1.  In  this  regime,  the  statistics  of  the 

displacement  and  inclination  of  the  shock  front  compare  favorably  with  the  linear 
analysis  predictions.  For  flows  with  Mi  >  M^'  —  1,  shock  waves  no  longer  have 
well-defined  fronts;  shock  wave  thickness  and  strength  vary  widely  in  the  transverse 
directions.  Multiple  peaks  in  pressure  are  found  along  the  mean  streamline  where 
the  local  thickness  of  the  shock  wave  has  increased  significantly. 
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CHAPTER  1 


INTRODUCTION 


1.1  Motivation 

A  fundamental  understanding  of  compressible  turbulence  is  necessary  for  the 
development  of  supersonic  transport  aircraft.  Compressibility  effects  on  turbulence 
were  found  significant  when  the  energy  associated  with  the  dilatational  fluctuations 
is  large  or  when  the  mean  flow  is  significantly  distorted —  expanded  or  compressed. 
The  presence  of  shock  waves  is  an  important  feature  that  distinguishes  high-speed 
flows  from  low-speed  ones.  Understanding  the  mechanisms  of  isotropic  turbulence 
interacting  with  a  shock  wave  is  not  only  of  generic  interest,  but  also  of  fundamental 
importance  in  understanding  the  interactions  of  turbulent  boundary  layers  with 
shock  waves  which  occur  in  many  practiced  engineering  applications;  the  flow  inside 
a  high  speed  compressor  or  a  gas  turbine,  the  flow  over  wings  in  supersonic  aircrafts, 
and  the  intake  flow  to  a  supersonic  ramjet  (scramjet). 

The  numericed  simulation  using  turbulence  models  is  becoming  a  standard  tool 
in  aerospace  technology.  Most  current  models  of  compressible  turbulence  are,  how¬ 
ever,  based  on  incompressible  turbulence  models.  A  better  understanding  of  the 
underlying  physics  could  lead  to  improvements  to  turbulence  models,  leading  to 
more  efficient  designs.  There  is,  therefore,  a  need  to  assess  our  understanding  of 
compressible  turbulence. 

The  present  work  is  a  fundamental  study  of  the  interactions  of  a  shock  wave  with 
turbulence.  We  investigate  the  interaction  of  isotropic  turbulence  with  a  shock  wave 
using  direct  numerical  simulation  and  linear  analyses. 


1.2  Survey  of  Previous  Work 

Studies  of  the  interaction  of  turbulence  with  a  shock  wave  were  initiated  using 
linear  theories  in  the  early  1950’s.  Twenty  years  later,  there  was  a  resurgence  of 
research  interest  in  this  area  through  experiments  on  tlie  interaction  of  isotropic 
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turbulence  and  turbulent  boundary  layer  with  a  shock  wave.  However,  direct  nu¬ 
merical  simulation  of  the  interaction  of  “true”  turbulence  with  a  shock  wave  has 
never  been  attempted. 

1.2.1  Linear  Analysis 

Using  his  general  theory  of  aerodynamic  sound  generation,  Lighthill  [1953]  esti¬ 
mated  the  acoustic  energy  scattered  from  the  interaction  of  turbulence  with  sound 
and  shock  waves.  Most  of  the  analytical  studies  of  shock-turbulence  interaction 
[Ribner  1953,  Moore  1953,  Kerrebrock  1956,  Chang  1957,  McKenzie,  and  Westphal 
1968]  are  based  on  the  linear  theories  of  three-dimensional  disturbances  interacting 
with  a  shock  wave.  These  disturbances  were  waves  of  vorticity,  entropy,  or  sound. 
Kovasznay[l953]  pointed  out  that  they  are  linearly  independent  in  weak  turbulence. 
Any  one  such  wave  interacting  with  the  shock  wave  generates  all  three  kinds  of  fluc¬ 
tuations  downstream  of  the  shock  wave.  The  linear  theories  developed  by  various 
researchers  followed  procedures  that  are  mathematically  different  but  physically 
equivalent  and  are,  therefore,  mutually  consistent:  inviscid  linear  equations  for  the 
disturbances  are  solved  downstream  of  the  shock,  and  the  boundary  conditions  at 
the  downstream  side  of  the  shock  front  are  expressed  in  terms  of  the  upstream 
disturbances  by  the  use  of  Rankine-Hugoniot  relations.  Ribner  [1953]  investigated 
the  passage  of  a  single  vorticity  wave  through  a  plane  shock  and  the  modification 
of  the  vortiiity  wave  with  simultaneous  generation  of  an  acoustically  intense  sound 
wave  in  a  reft  runce  frame  fixed  on  the  shock  wave.  He  later  extended  this  analysis 
to  study  turbulence  amplification  due  to  a  shock  wave  [1954]  and  the  flux  of  acous¬ 
tic  energy  emanating  on  the  downstream  side  of  the  shock  [1969].  He  also  used  it 
to  predict  the  one-dimensional  power  spectra  of  various  fluctuations  downstream 
of  the  shock  [1987].  Moore  [1953]  analyzed  the  flow  field  produced  by  the  obhque 
impingement  of  weak  plane  disturbances  on  a  normal  shock  wave  in  a  reference 
frame  fixed  on  the  mean  upstream  flow.  Chang  [1955]  investigated  the  interaction 
of  a  plane  shock  and  oblique  plane  flisturbances  with  special  reference  to  entropy 
w-aves.  McKenzie  and  VVestph;d  [1968]  investigated  the  effect  of  a  stationary  plane 
shock  on  the  travelling  waves  of  vorticity,  entropy,  and  sound.  Anyiwo  and  Bushnell 
'1982]  revisited  the  analysis  of  McKenzie  and  Westphal  [1968]  to  identify  primary 
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mechanisms  of  turbulence  enhancement —  amplification  of  vorticity  mode,  genera¬ 
tion  of  vorticity  mode  from  the  interaction  of  acoustic  and  entropy  modes  with  a 
shock  wave,  and  turbulence  “pumping”  by  shock  oscillations. 

Debieve,  Gouin,  and  Gaviglio  [1982a,  1982b]  analyzed  turbulence  evolution  through 
the  shock  using  the  Reynolds  stress  transport  equation.  They  were  able  to  separate 
the  effects  of  the  specific  turbulent  sources  from  the  effects  of  the  mean  motion — 
convection  and  oroduction.  Their  prediction  of  the  longitudinal  velocity  fluctuation 
showed  good  comparison  with  the  experimental  result. 

1.2.2  Experiments 

There  has  been  a  significant  accumulation  of  experimental  data  on  the  shock  tur¬ 
bulence  interaction  during  the  last  decade.  Debieve,  Gouin,  and  Gaviglio  [1982a, 
1982b|  performed  an  experiment  on  the  turbulent  boundary  layer  interacting  with  a 
shock  wave.  They  measured  the  mean  and  turbulent  fields  in  an  adiabatic  compres¬ 
sion  ramp,  where  the  mean  upstream  Mach  number  was  =  2.32,  with  a  corner 

angle  of  6°.  They  found  amplifications  of  turbulence  intensity  Uj  ,  the  structure 

parameter  -u\u2/u^  ,  and  the  temperature  fluctuation  inside  the  boundary  layer. 

Dolling  and  Or  [1985]  measured  wall  pressure  fluctuations  upstream  of  the  corner 
in  flows  with  =  3.0  over  compression  ramps  with  corner  angles  of  8°,  12°,  16°, 
and  20°.  They  found  that  the  shock  wave  structure  is  unsteady  in  both  separated 
and  attached  downstream  flows,  resulting  in  a  region  in  which  the  wall  pressure 
signal  is  intermittent. 

Andreopoulos  and  Muck  [1987]  investigated  the  wall  pressure  fluctuations  in  the 
shock-wave/boundary-layer  interactions  over  two-dimensional  ramps,  and  found 
that  the  frequency  of  the  shock-wave  unsteadiness  is  of  the  same  order  as  the  burst¬ 
ing  frequency  of  the  upstream  boundary  layer  and  independent  of  the  downstream 
separated  flow. 

Smits  and  Muck  [1987]  performed  experiments  to  study  the  effects  of  different 
compression  corners  with  angles  of  8°,  16°,  and  20°  on  a  compressible  turbulent 
boundary  layer  with  =  2.9.  They  found  that  the  interaction  significantly  am¬ 
plifies  the  turbulent  stresses,  and  that  the  amplification  increases  with  increasing 
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the  effective  normal  Mach  number.  The  structure  parameter  —u'^u'2/ also  in¬ 
creased  significantly,  which  was  attributed  mainly  to  the  unsteady  oscillation  of  the 
shock  system. 

Kuntz,  Amatucci,  and  Addy  [1987]  conducted  an  experimental  investigation  of 
the  interaction  between  a  shock  wave  and  a  turbulent  boundary  layer.  Compression 
corners  were  used  to  generate  an  oblique  shock  wave  in  the  flow  field  with  = 
2.94.  Ramp  angles  of  8°,  12°,  16°, 20°,  and  24°  were  used  to  produce  a  range 
of  possible  flow  fields,  including  flows  with  no  separation,  incipient  separation, 
and  significant  separation.  They  found  that  that  the  boundary  layer  after  the 
interactions  showed  an  acceleration  of  the  mean  flow  near  the  wall  as  the  boundary 
layers  began  to  return  to  equilibrium,  and  that  the  mean  streamwise  velocity  profiles 
downstream  of  the  separated  compression  corner  were  wavy  due  to  the  redeveloping 
boundary  layer  which  had  a  velocity  profile  with  inflection  at  reattachment. 

There  are  a  number  of  additional  experimental  studies  of  the  interaction  of  tur¬ 
bulent  boundary  layers  with  an  oblique  shock  [Settles,  Fitzpatrick,  and  Bogdonoff 
1979,  Dussauge,  Muck,  and  Andreopoulos  1986,  Jayaram,  Taylor,  and  Smits  1987, 
Selig,  Andreopoulos,  Muck,  Dussauge,  and  Smits  1989].  A  general  finding  from 
these  experiments  is  that  Reynolds  shear  stress  and  turbulence  intensities  are  am¬ 
plified  across  the  shock  wave.  The  studies  of  oblique  shock  wave/turbulent  bound¬ 
ary  layer  interaction  included  several  additional  phenomena  which  complicated  the 
flow  behavior.  These  phenomena  are:  (a)  oscillation  of  the  shock  wave  in  the  lon¬ 
gitudinal  direction,  (b)  flow  separation  downstream  of  the  shock,  (c)  streamline 
curvature,  and  (d)  wall  effects  which  result  in  high  turbulence  intensity  and  high 
flow  anisotropy.  Because  of  these  complications,  it  was  impossible  to  identify  the 
sole  effect  of  a  shock  wave  on  turbulence. 

In  order  to  isolate  the  effects  of  a  shock  wave  on  turbulence,  several  experiments 
on  the  interaction  between  the  shock  wave  and  grid-generated  turbulence  have 
been  performed.  Debieve  and  Lacharme  [1986]  experimentally  investigated  a  shock- 
wave/free  turbulence  interaction  at  ~  2.3  over  a  ramp  with  a  corner  angle  of 
6°.  They  njeasured  velocity  and  temperature  spectra  upstream  and  downstream 
of  the  shock  wave  and  concluded  that  turbulent  fluctuations  are  amplified  and  the 
Taylor  microscales  increase  during  the  interaction.  An  intermittency  effect  due  to 
unsteady  shock  wave  distortion  on  turbulence  statistics  was  also  clearly  described. 
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Keller  and  Merzkirch  [1990]  performed  an  experiment  on  the  interaction  of  grid 
generated  turbulence  with  a  shock  inside  a  shock  tube.  They  verified  amplification 
of  the  turbulence  intensity  quantitatively,  showing  that  amplification  was  restricted 
to  the  lower  wave  numbers  in  the  spectrum.  This  was  consistent  with  the  conclusion 
of  length  sccde  increase  made  by  Debieve  and  Lacherme  [1986],  but  it  contradicts  the 
intuitive  reasoning  that  mean  flow  compression  decreases  the  relevant  turbulence 
length  scedes.  Our  results  (Sec.  2.1  and  Subsection  4.1.4)  show  that  the  Taylor 
microscale  actually  decreases  in  passing  through  the  shock. 

Honkan  and  Andreopoulos  [1990]  examined  the  interaction  of  a  normal  shock 
wave  with  homogeneous  grid-generated  turbulence.  They  found  that  turbulence 
is  considerably  amplified  during  the  interaction,  and  that  the  amplification  ratio 
of  turbulence  is  not  the  same  for  different  length  scales  and  different  turbulence 
intensities. 

Jacquin,  Blin  and  Geffroy  [1991]  investigated  the  interactions  of  a  normal  shock 
wave  with  grid-generated  turbulence  and  a  turbulent  jet,  ?*nd  compared  turbulence 
amplifications  with  the  predictions  by  a  linear  analysis.  They  observed  that  turbu¬ 
lence  amplification  was  not  significant  during  the  interaction,  and  that  the  decay 
of  turbulent  kinetic  energy  was  accelerated  downstream  of  the  shock  wave. 

The  aforementioned  experiments  treated  the  interaction  of  a  shock  with  quasi- 
incompressible  turbulence  where  fluctuations  in  pressure  and  density  are  not  sig¬ 
nificant.  A  comprehensive  experiment  on  the  interaction  of  weak  shocks  {M^  — 
1.007,1.03,  and  1.1)  with  a  random  medium  of  density  inhomogeneity  was  per¬ 
formed  by  Hesselink  and  Sturtevant  [1988].  They  observed  that  the  pressure  his¬ 
tories  of  the  distorted  shock  waves  were  both  peaked  and  rounded.  In  the  rounded 
case,  they  found  the  perturbed  shock  was  made  up  of  a  succession  of  weak,  slightly 
curved  fronts,  and  the  total  effective  shock  thickness  was  significantly  greater  than 
the  classical  Taylor  thickness.  They  concluded  that  the  observed  distortions  of  the 
shock  can  best  be  explained  in  terms  of  the  focusing/defocusing  of  its  front  due  to 
inhomogeneity  of  the  medium. 

1.2.3  Numerical  Simulation 

Zang  et  al.  [1984]  simulated  the  interaction  of  a  shock  wave  with  a  single  wave 
using  two-dimensional  Euler  equations.  The  interaction  of  a  shock  wave  with  an 
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isolated  flow  inhomogeneity  was  computed  by  Hussaini  et  al.  [1986]  and  Meadows 
et  al.  [1991],  and  with  two-dimensional  turbulence  by  Rotman  [1991]  and  Lee,  Lele, 
and  Moin  [1991a]. 

Rotman  [1991]  numerically  calculated  the  change  in  a  turbulent  flow  caused  by 
the  passage  of  a  travelling  shock  wave.  He  found  that  the  shock  causes  an  increase  in 
the  turbulent  kinetic  energy  and  that  the  length  scale  of  the  turbulent  field  behind 
the  shock  is  smaller  than  that  in  front.  He  also  found  that  increasing  the  initial 
turbulent  kinetic  energy  caused  a  straight  shock  wave  to  evolve  into  a  distorted 
front. 

Lee,  Lele,  and  Moin  [1991a]  found  that  vorticity  amplification  in  the  numerical 
simulation  compared  well  with  the  predictions  of  the  linear  analyses,  but  that 
turbulent  kinetic  energy  evolution  behind  the  shock  showed  significant  nonlinear 
effects.  The  energy  spectrum  was  found  to  be  enhanced  more  at  small  scales, 
leading  to  an  overall  length  scale  decrease. 

Zang  et  al.  [1984]  examined  various  effects  pertinent  to  the  amplification  and 
generation  of  turbulence  in  shock/turbulent  boundary  layer  interaction  and  placed 
limits  on  the  range  of  validity  of  linear  theory.  Hussaini  et  al.  [1986]  numerically 
investigated  the  effects  of  upstream  eddy  motion  and  temperature  inhomogeneity  on 
the  enhancement  and  production  of  turbulence.  Meadows  et  al.  [1991]  computed 
two-dimensional  shock-vortex  interaction  using  a  shock  capturing  scheme.  They 
qualitatively  evaluated  the  effects  of  upstream  vortex  strength  on  both  the  flow 
field  and  acoustic  field  generated  by  the  interaction. 


1.3  Objectives  and  Overview 

The  primary  objective  of  this  work  is  to  investigate  the  physics  of  the  interaction 
of  isotropic  turbulence  with  shock  waves  using  direct  numerical  simulations  and 
linear  analyses.  The  simulations  and  linear  analyses  provide  statistical  information 
for  testing  of  turbulence  models.  Instantaneous  flow  fields  from  the  simulations 
contribute  to  our  understanding  of  the  physical  nature  of  the  interaction. 

The  principal  contributions  and  findings  of  this  work  are  as  follows: 
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•  A  numerical  scheme  to  generate  turbulence  at  the  inflow  boundary  was  devel¬ 
oped  and  the  simulation  of  spatially  evolving  grid-generated  turbulence  was  con¬ 
ducted. 

•  The  linear  mechanisms  involved  in  the  shock/turbulence  interaction  enhance  the 
amplitude  of  turbulent  fluctuations  and  decrease  the  relevant  turbulence  length 
scales. 

•  Linear  ancdysis  predicts  the  corrugation  of  the  shock  front  caused  by  upstream 
turbulence  fluctuations.  The  distortion  of  the  shock  front  is  scaled  with  the 
upstream  turbulence  intensity  and  length  scales. 

•  Isotropic  upstream  turbulence  becomes  axisymmetric  after  the  interaction.  All 
the  components  of  turbulent  kinetic  energy  are  amplified  across  the  shock  wave. 
The  streamwise  intensity  is  amplified  more  than  the  transverse  components. 
Fluctuations  in  pressure,  density,  and  temperature  are  significantly  enhanced. 

•  Power  spectra  of  turbulent  fluctuations  are  more  amplified  at  small  scales  than 
at  large  scales.  The  integral  turbulence  length  scale  and  Taylor  microscales 
decrease  during  the  interaction. 

•  Transverse  components  of  vorticity  are  enhanced  because  of  the  mean  flow  com¬ 
pression,  but  the  component  normal  to  the  shock  remains  unchanged.  Baro- 
clinic  torque  has  a  negligible  contribution  to  the  production  of  vorticity  in  the 
shock/quasi-incompressible  turbulence  interaction. 

•  Rapid  evolution  of  turbulent  kinetic  energy  found  downstream  of  the  shock  wave 
is  caused  by  the  nonlinear  pressure  work.  Decomposition  of  the  pressure  work 
term  shows  that  the  inhomogeneous  pressure  transport  is  the  main  cause  of  the 
rapid  evolution. 

•  Isentropic  relations  hold  between  normalized  fluctuations  in  pressure,  density, 
and  temperature  throughout  the  flow  field. 

•  Shock  waves  are  found  to  be  distorted  by  the  upstream  turbulence.  Instanta¬ 
neous  shock  wave  structure  depends  on  the  upstream  fluctuation  Mach  number 
and  the  mean  shock  strength,  k^or  flows  with  small  fluctuation  Mach  numbers, 
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shock  waves  have  a  well-defined  front.  In  this  regime,  the  statistics  of  the  dis¬ 
torted  shock  front  compare  favorably  with  the  linear  analysis  predictions.  For 
flows  with  large  fluctuation  Mach  numbers,  multiple  peaks  in  pressure  are  found 
along  the  mean  streamline  where  the  local  thickness  of  a  shock  wave  has  increased 
significantly. 

This  report  is  organized  as  follows.  In  Chapter  2  linear  analyses  of  shock-isotropic 
turbulence  interaction  are  conducted.  The  governing  equations  and  the  numerical 
method  chosen  are  discussed  in  Chapter  3.  Chapter  4  describes  the  numerical 
simulations  of  shock-turbulence  interaction  and  their  results.  Conclusions  and  rec¬ 
ommendations  for  future  work  are  given  in  Chapter  5.  The  appendices  include  brief 
descriptions  of  the  linear  analyses  (Appendix  A  and  B),  a  proposal  for  an  alias-free 
compressible  turbulence  simulation  (Appendix  C),  development  and  validation  of 
a  numerical  method  for  the  simulations  of  spatially  evolving  turbulence  (Appendix 
D),  the  limitations  on  the  physical  parameters  for  direct  numerical  simulation  of 
shock-turbulence  interaction  (Appendix  E),  the  effect  of  the  shock  oscillation  on 
turbulence  statistics  (Appendix  F),  the  effect  of  outflow  boundary  conditions  (Ap¬ 
pendix  G),  and  the  drift  in  the  mean  shock  position  and  the  outflow  condition 
(Appendix  H). 
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CHAPTER  2 


LINEAR  ANALYSIS 

Two  different  linear  approaches  are  used  to  investigate  the  interaction  of  isotropic 
turbulence  with  a  normal  shock  wave.  The  first  approach  is  the  rapid  distortion 
theory  (RDT).  Th^  second  is  the  compressible  linear  interaction  analysis  (LIA). 
Both  analyses  are  for  inviscid  flows.  The  assumptions  and  main  features  of  these 
analyses  are  discussed  in  Appendices  A  and  B.  The  main  effect  of  a  shock  wave  on 
turbulence  is  the  mean  flow  compression  in  the  direction  normal  to  the  shock  wave. 
The  secondary  effects  are  the  vorticity  generation  due  to  shock  front  curvature 
and  turbulent  kinetic  energy  generation  caused  by  the  unsteady  movement  of  the 
shock  front.  RDT  accounts  for  the  effect  of  the  mean  flow  compression,  while  LIA 
includes  all  three  effects. 

For  proper  applications  of  the  linear  analyses,  upstream  Mach  number  variation 
may  be  considered  as  a  small  perturbation  from  the  mean  Mach  number.  Fur¬ 
thermore,  time  required  for  turbulence  to  pass  throu  ;h  the  shock  wave  may  be 
considered  small  compared  to  a  turbulence  time  scale,  so  that  turbulence  has  no 
time  to  redistribute  energy  into  different  scales  through  nonlinear  processes. 

In  the  following,  interaction  of  isotropic  turbulence  with  a  normal  shock  wave  is 
discussed.  Because  LIA  is  more  comprehensive,  its  results  are  mentioned  in  detail, 
and  those  from  RDT  are  introduced  for  comparison  when  needed.  Turbulence  be¬ 
fore  the  interaction  is  considered  to  be  purely  vortical,  that  is,  upstream  turbulence 
has  no  fluctuations  in  density  and  entropy.  The  fluid  is  assumed  to  be  ideal  gas 
with  the  specific  heat  ratio  7  =  1.40. 

The  coordinate  system  used  in  the  analysis  is  Cartesian,  as  shown  in  Figure 
2.1.  The  streamwise  direction  xj  is  aligned  with  th^  direction  normal  to  the  mean 
shock  plane.  The  reference  frame  is  fixed  on  the  mean  shock  plane.  In  this  system, 
upstream  flow  approaches  the  shock  with  a  supersonic  speed,  and  downstream  flow 
leaves  with  a  subsonic  speed. 


2.1  Turbulence  Statistics 


Figure  2.2  presents  amplifications  of  the  transverse  vorticity  components  pre¬ 
dicted  by  LIA  and  RDT.  The  streamwise  vorticity  component  is  unchanged  through 
the  linear  interaction.  As  the  strength  of  the  shock  wave  increases  (higher  mean 
Mach  number),  the  ratio  of  downstream  to  upstream  vorticity  also  increases.  We 
find  that  the  asymptotic  value  of  the  amplification  factor  for  mean  square  vorticity 
for  a  shock  wave  with  very  large  Mach  number  is  about  20  in  LIA  and  36  in  RDT. 
The  ratio  of  the  vorticity  amplification  by  RDT  is  simply  the  ratio  of  downstream 
to  upstream  density,  which  is  explained  as  an  enhancement  of  vorticity  due  to  the 
shrinking  of  the  cross  section  of  a  transversely-oriented  vortex  tube  by  the  mean 
flow  compression.  The  predictions  by  LIA  and  RDT  agree  very  well  for  weak  shock 
waves  but  do  not  compare  as  well  for  stronger  shock  waves.  This  implies  that  ef¬ 
fects  of  shock  front  curvature  and  shock  front  unsteadiness  are  negligible  for  weak 
shock  waves,  while  these  secondary  effects  become  more  significant  for  stronger 
shock  waves.  Prediction  of  lower  amplification  by  LIA  suggests  that  secondary 
mechanisms  have  adverse  effects  on  enhancement  of  vorticity  fluctuations. 

Figure  2.3  shows  the  amplification  of  solenoidal  turbulent  kinetic  energy  by  LIA 
and  RDT.  Both  approaches  predict  more  enhanced  streamwise  fluctuations  than 
spanwise  fluctuations  for  shock  waves  with  the  mean  upstream  Mach  number  < 
2.0.  However,  this  trend  is  reversed  for  stronger  shock  waves  in  LIA  predictions.  In 
fact,  the  predictions  by  LIA  and  RDT  are  close  only  for  very  weak  shocks,  differing 
significantly  for  =-  1.5. 

Through  interaction  of  vortical  waves  with  a  shock  wave,  acoustic  waves  are  gen¬ 
erated  downstream  of  the  shock  wave.  These  acoustic  waves  accompany  the  purely 
dilatational  velocity  fluctuations,  which  also  contribute  to  the  total  turbulent  ki¬ 
netic  energy.  LIA  can  also  predict  the  acoustic  energy  generated  downstream  of 
the  shock.  Figure  2.4  presents  the  amplification  of  turbulent  kinetic  energy,  in¬ 
cluding  both  solenoidal  (or  vortical,  incompressible)  and  dilatational  (or  acoustic, 
compressible)  velocity  fluctuations.  Since  part  of  the  acoustic  energy  undergoes 
an  inviscid  decay  (ref.  Appendix  B),  turbulence  behind  the  shock  is  not  homoge¬ 
neous  in  the  streamwise  direction.  Figure  2.4  presents  velocity  fluctuation  levels 
at  both  immediate  downstream  (near-field)  and  far  downstream  (far-field)  of  the 
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shock.  Streamwise  velocity  fluctuation  is  always  larger  than  that  of  the  spanwise 
velocity  fluctuations  in  the  near-field.  The  far-field  velocity  fluctuations  are  com¬ 
posed  mostly  of  the  solenoidal  velocity  fluctuations:  the  acoustic  waves  contribute 
less  than  2%  of  the  far-field  turbulent  kinetic  energy  for  =  1.2. 

According  to  LI  A,  the  vorticity  waves  incident  at  angles  beyond  a  critical  angle 
&cr  =  )  generate  acoustic  waves  which  decay  exponentially  as  they  prop¬ 

agate  downstream.  This  leads  to  an  inviscid  decay  of  the  compressible  part  of 
velocity  fluctuations.  The  decays  of  velocity  fluctuations  are  shown  in  Figure  2.5 
for  upstream  turbulence  with  a  spectrum  of 

f:(*)~(^)‘‘exp|-2(^)2|,  (2.1) 

iCo  Kq 

where  ko  is  the  characteristic  wave  number  corresponding  to  the  energy  peak.  This 
is  the  form  of  the  inflow  turbulence  spectrum  used  in  the  direct  numerical  simulation 
of  shock  turbulence  interaction  in  Chapter  3.  Significant  but  monotonic  decay 
occurs  just  downstream  of  a  shock  wave  which  is  caused  by  the  inviscid  decay  of 
acoustic  waves. 

Decays  of  velocity  fluctuations  for  different  shock  strengths  are  shown  in  Figure 
2.6.  The  decay  in  the  downstream  velocity  fluctuation  is  monotonic  for  all  shock 
strengths. 

We  also  investigated  the  effect  of  upstream  spectrum  shape  on  the  decay,  as 
shown  in  Figure  2.7.  In  addition  to  the  spectrum  described  in  (2.1),  we  used  the 
von  Karman  spectrum 


E{k) 


{kfkoY 


(2.2) 


This  is  a  good  approximation  for  the  high  Reynolds  number  turbulence,  where  the 
slope  is  k'^  at  small  k  and  k~^^^  for  the  inertial  subrange  at  large  k.  Monotonic 
decay  of  turbulent  kinetic  energy  is  reproduced.  However,  the  decay  rate  is  faster 
compared  to  that  of  (2.1)  (see  Figure  2.5),  since  the  von  Karman  spectrum  has 
more  small  scale  content,  thus  causing  faster  decay. 
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Experimental  studies  [Debieve  et  al.  1986,  Keller  et  al.  1990]  have  reported 
that  large  scale  turbulent  motions  are  enhanced  more  than  small  scale  motions  as 
turbulence  passes  through  a  shock  wave,  leading  to  the  overall  increase  of  turbulence 
length  scales.  In  order  to  check  whether  turbulence  length  scales  do  indeed  increase, 
we  investigate  the  amplification  of  the  one-dimensional  velocity  spectrum  Ei(fci), 
which  is  defined  as 


/  +  00  ^-foo 

/  En(k)dfc2^fc3,  (2.3) 

-oo  J  —oo 

where  k  =  is  the  wave  number  vector.  The  velocity  spectrum  tensor 

Eij(k}  is  defined  as 


(2.4a) 

where  (•)  denotes  the  Fourier  transform,  the  superscript  *  denotes  complex  conju¬ 
gate,  and  (•)  indicates  the  ensemble  average.  For  incompressible  isotropic  turbu¬ 
lence 


^ij  ( k  )  —  > 


where  k  =  |k|,  and  E(k)  is  the  energy  spectrum  function. 

In  the  following  analysis,  we  consider  only  contributions  from  the  vorticity  waves 
to  the  upstream  and  downstream  velocity  fluctuations.  A  velocity  fluctuation  u 
associated  with  a  vorticity  wave  in  the  homogeneous  field  can  be  represented  as 

u  =  u(u;,^'2,A:3;  Ei)exp[i(k  •  x  -  u;/)],  (2.5) 

where  k  •  u  =  0  for  a  solenoidal  wave.  Assuming  that  vorticity  fluctuations  are 
simply  advected  by  the  mean  flow,  we  obtain  the  following  dispersion  relation 

u;-f/,iti=0.  (2.6) 

In  this  approximation,  the  turbulence  fields  upstream  and  downstream  of  the  shock 
are  considered  homogeneous  and  frozen  with  respect  to  the  corresponding  mean 
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flow.  From  Appendices  A  and  B,  the  Fourier  coefficient  of  the  streamwise  veloc¬ 
ity  fluctuation  after  the  interaction  can  be  expressed  in  terms  of  that  before  the 
interaction.  The  transfer  function  Z(a;,fc2,A:3)  is  defined  as 


Z{u;,k2,k3;M^) 


(2.7) 


where  the  superscripts  U  and  D  refer  to  upstream  and  downstream  states. 

If  there  is  no  generation  or  destruction  of  waves  inside  the  shock  wave,  the 
frequency  u;  of  a  vorticity  wave  remains  unchanged  by  the  interaction,  while  the 
associated  wave  number  ki  changes  to  satisfy  the  dispersion  relation  (2.6).  The 
one-dimensional  frequency  spectrum  E\{u})  is  defined  as 


Ei{uj)  =  J  J  Fii(a},k2,k3;Ui)dk2dk2, 


(2.8a) 


where  the  limits  of  the  integrations  are  (  — oo,  -f-oo)  for  both  ^2  and  ^-3.  In  the  above 
expression,  Fn  is 


Fii{ui,k2,ky,Ui)  =  ui{u;,k2yk:i-,Ui)u’^{u),k2,k:i;Ui).  (2.86) 

Amplification  of  the  one-dimensional  frequency  spectrum  is  described  by  the 
ratio  of  the  frequency  spectra  before  and  after  the  interaction  5‘*'(u»;  AI^ ),  which  is 


5‘"(u;;Mf) 


J  J  F{iiu;,k2,kr,Uy)dk2dk:i 
JjFl[{u>,k2,kr,U[^)dk2dk^ 


J  J\Z{u,yk2,k3;M!;^)\^Fu{^,k2,kr,U{^)dk2dk2 

1 1  En{^,h,h;UP)dk2dk3 


(2.9) 


Spectral  amplification  ratio  is  dependent  on  the  special  form  of  the  upstream  spec¬ 
trum  shapes.  We  choose  a  von  Karnian  spectrum  (2.2)  for  the  energy  spectrum  in 
the  following  analyses. 


The  amplification  ratios  of  the  one-dimensional  frequency  spectra  for  different 
shock  strengths  are  shown  in  Figures  2.8:  the  results  from  RDT  and  LIA  are  pre¬ 
sented  in  Figures  2.8(a)  and  2.8(b),  respectively.  The  results  from  RDT  and  LIA 
are  qualitatively  the  same:  the  spectrum  amplification  ratio  is  larger  for  a  wave 
with  small  u;,  which  is  consistent  with  the  results  of  Ribner  [1987].  Some  researchers 
interpret  this  fact  as  evidence  of  turbulence  length  scale  increase  through  the  inter¬ 
action,  but  this  conclusion  does  not  necessarily  follow,  for  the  change  in  frequency 
spectra  reflects  a  change  in  time  scale,  not  in  length  scale.  More  amplification  at  the 
small  frequency  part  of  the  spectrum  implies  that  turbulence  time  scale  increases 
through  the  interaction. 

To  investigate  the  length  scale  change,  one  should  evaluate  the  amplification  of 
wave  number  spectrum  rather  than  that  of  frequency  spectrum.  Amphfications  of 
these  two  spectra  are  not  the  same:  as  a  wave  passes  through  the  shock,  the  relation 
between  wave  number  and  frequency  changes  due  to  the  mean  flow  deceleration  (see 
(2.6)).  The  one-dimensional  wave  number  spectrum  is  obtained  by  replacing  uj  with 
Uiki  in  (2.8b).  Streamwise  velocity  fluctuations  upstream  and  downstream  of  the 
shock  can  be  represented  as 

uf  =  J[J  J  Fl^i{Ul'ki,k2,h;U\')dk2dh  dk^,  (2.10a) 

/  /  Fn{uFki,k2,ky,uP)dk2dk:^  dk^ .  (2.106) 

Downstream  velocity  fluctuation  can  also  be  expressed  in  terms  of  the  upstream 
spectrum  and  its  amplification  ratio  at  the  same  wave  number  across  the  shock 
wave  5^*(fci ;  )  as 

j  S^^(kyM';)^j  j  F\[(U\’kukykyU\^)dk2dk^  dfcf 
S^^(kyM\)^j  j  F\\(V^  kyk2,kyUhdk2dk2  (Jdk^), 

(2.11a) 
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where 


_  s s F{>(U{>ki,k2,ky,U{’)‘‘h<th 
'  SjFj>^{V’;ki,k2,k3-,U[')dk2dk3' 

Here  J  is  the  Jacobian  of  the  transformation  from  the  upstream  to  the  downstream 
wave  number  definea  as  J  =  <  1.  (The  upstream  wave  number  interval 

dk^  corresponds  to  the  downstream  wave  number  interval  dk^  =  J~^dk^ .)  Com¬ 
paring  the  terms  in  (2.106)  and  (2.11a),  we  have  the  amplification  ratio  of  the 
one-dimensional  wave  number  spectrum  JS^^{ki;  )  as 


=  J  X 


J  J  Ff;(Uf>k,,k2,k,;  Uf)dk2dk3 

f  f  Ff'(  JUf  ».-2,  tj;  (/''  )dk2dk3 


S  S  F"{JU"  k,,k2,k3-.Ui')dk2dk3 
S  S  F!^,{U« k3,k2,k3-,U’^)dk2dk3 


=  J  X  S"(a.;Mf')  X 


E^(u,/J)’ 


(2.12) 


where  u  =  UPki  =  ju¥ ki-  The  wave  number  spectrum  amplification  factor  is, 
therefore,  the  product  of  the  Jacobian,  the  frequency  spectrum  amplification  and 
the  ratio  of  the  upstream  frequency  spectra  at  two  different  frequencies.  The  resul¬ 
tant  wave  number  spectrum  amplifications  for  different  shock  strengths  predicted 
by  RDT  and  LIA  are  presented  in  Figures  2, 9(a)  and  (b),  respectively.  These  pre¬ 
dictions  are  qualitatively  consistent:  the  spectral  density  increases  more  at  large 
wave  numbers,  even  though  it  increases  more  at  small  frequencies  (see  Figure  2.8). 
LIA  predicts  a  suppression  of  the  spectrum  at  small  wave  numbers  for  large  up¬ 
stream  Mach  numbers,  while  RDT  predicts  an  amplification  at  all  wave  numbers 
irrespective  of  the  shock  strength.  Larger  amplification  at  large  wave  numbers  is 
more  pronounced  for  stronger  shock  waves. 

It  is,  therefore,  erroneous  to  infer  an  increase  in  the  length  scale  in  shock  turbu¬ 
lence  interaction  by  appealing  to  Ribner’s  analysis  (akin  to  Figures  2.8)  as  Keller 
et  al.  [1990]  have  done.  Investigation  of  the  spectrum  amplification  leads  to  the 
conclusion  that  turbulence  time  scale  increases  through  the  interaction,  while  tur¬ 
bulence  length  scale  decreases. 
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The  experimental  results  by  Debieve  et  al.  [1986]  are  consistent  with  the  present 
predictions.  However,  they  compared  the  upstream  frequency  spectrum  to  that 
on  the  shock  and  concluded  an  increase  in  the  length  scale.  The  spectrum  on 
the  shock  is  contaminated  by  the  intermittency  effect  due  to  the  unsteady  shock 
front  distortion.  The  characteristic  length  scale  of  the  distortion  is  scaled  with  the 
upstream  turbulence  length  scale  (see  Section  2.2  for  LIA  and  Section  4.2  for  DNS) 
to  yield  an  apparent  amplification  of  the  spectrum  on  the  shock  at  an  energetic 
(or  large)  scale.  This  enhancement  does  not  necessarily  imply  the  amplification 
of  the  turbulent  motion  of  that  scale.  In  order  to  investigate  length  scale  change, 
one  has  to  transform  frequency  spectra  into  wave  number  spectra  and  compare 
the  upstream  spectrum  with  the  downstream  spectrum.  Proper  comparison  of 
wave  number  spectra  shows  more  amplification  at  small  scales  rather  than  at  large 
scales,  leading  to  an  overall  scale  decrease.  Direct  numerical  simulations  confirm 
that  turbulence  length  scales  do  decrease  through  the  shock-turbulence  interaction 
(see  Section  4,2). 

Since  linear  analyses  predict  the  change  of  the  Fourier  coefficients  of  the  velocity 
components  across  the  shock  wave  (see  (2.7)),  the  changes  in  turbulence  length 
scales  can  also  be  calculated.  The  change  in  the  Taylor  microscale  Aq  defined  as 


(2,13) 

across  the  shock  wave  is  independent  of  the  shape  of  the  three-dimensional  energy 
spectrum  for  isotropic  upstream  turbulence,  since  the  contributions  to  the  integral 
involving  the  w’ave  number  magnitude  k  (see  (/1. 2.5))  are  cancelled  out  in  the  eval¬ 
uation  of  the  ratios  of  and  ^  across  the  shock  wave.  The  denominator  in 
(2.1.3),  can  be  evaluated  usiiie; 

/// (2.14) 

P’iKures  2.10(a)  and  (b)  show  the  change  in  the  Taylor  microscales  across  the  shock 
wave  predicted  by  LIA  and  HDT,  (In  Ll,\  prediction,  only  the  solenoidal  velocity 
compionent  is  included  because  the  contrilmtion  from  the  dilatational  component 


>^i-- 


u 
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becomes  negligible  after  a  short  distance  from  the  shock  wave.  See  Figures  2.4, 
2.5  and  2.6.)  All  the  length  scales  are  found  to  decrease  through  the  interaction: 
the  streamwise  sccde  decreases  more  than  the  transverse  scale.  LIA  predicts  more 
reduction  in  the  streamwise  scale  and  less  reduction  in  the  transverse  scale  than 
RDT  does.  For  weak  shock  waves,  however,  the  predictions  by  LIA  and  RDT  are 
in  good  agreement. 

Shock/turbulence  interaction  leads  to  noise  generation  behind  the  shock  wave  in 
the  form  of  fluctuating  pressure,  p^  Using  (B.55)  and  (R.56)  and  from  the  isotropy 
relations  (2.4)  for  upstream  velocity  fluctuations  (see  also  (j5.61)  -  (5.64)), 


ll  =  ±  f  V  !  ^ 

pD^  \  (7  +  1  )m  —  (7  —  1 )  /  J  cos^  0  cos^  ^ 


1 


(h  +  i) 


27771 

771  -  (7  -  1) 


700  7277 

/  Eik)dk 
Jo  Jo 


cos  9 

cos^  O' 


dO 


3  ^  / _ 27771 _ \  ^ 

8  Uf  \  (7  +  l)'^  -  (7  -  1)/  Jo  COS-0' 


(2.15) 


where  is  the  downstream  mean  pressure  and 

=  tan~‘(77i  tan  ^).  {B.5c) 

Figures  2.11(a)  and  (b)  show  the  pressure  fluctuation  for  various  shock  strengths 
both  at  the  immediate  downstream  and  far  downstream  of  the  shock  wave  normal¬ 
ized  by  the  downstream  mean  pressure  and  upstream  mean  pressure,  respectively. 
Note  that  for  normalization  we  have  also  used  the  upstream  turbulence  intensity, 
because  pressure  fluctuations  scale  with  the  upstream  turbulence  mtensity  inde¬ 
pendent  of  the  shape  of  the  spectrum.  Near  field  noise  scaled  with  the  downstream 
mean  pressure  peaks  around  A/j^  =  1.3,  and  far  field  noise  reaches  its  asymptoti¬ 
cally  maximum  strength  for  very  large  Mach  numbers.  Note  that  acoustic  energy 
decays  by  an  order  of  magnitude  from  the  near  field  to  the  far  field.  Even  though 
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the  near-  and  far-field  pressure  intensities  scaled  with  the  downstream  mean  pres¬ 
sure  tend  to  asymptotic  values  for  an  infinite  strength  shock  wave,  their  absolute 
amplitudes  increase  indefinitely  as  the  shock  strength  increases  (see  Figure  2.11(b)). 

2.2  Statistics  of  the  Shock  Front 

Using  LIA,  we  can  estimate  the  level  of  fluctuations  of  the  shock  front  caused 
by  the  action  of  turbulence.  The  details  for  calculating  the  variances  of  a  local 
shock  front  displacement  its  inchnation  angle  <72,  and  its  curvature  k,2  are  shown 
in  Appendix  B.  The  nondimensionalized  variances  of  those  quantities  arc 
and  where  ko  is  the  wave  number  corresponding  to  the  energy  peak  in  the 

spectrum. 

The  dimensionless  variance  of  the  shock  front  displacement  (see  (B.68))  can  be 
expressed  as 


ko.. 
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E{k)dk 
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-f  6^)  cosOdO, 


(2.16) 


where  k*  =  kjko,  and  Up  is  the  rms  fluctuation  velocity  in  one  direction.  (Defini¬ 
tions  of  ag  and  hg  are  given  in  (BAS)  and  (F.21).)  Likewise,  the  dimensionless 
variance  of  the  shock  front  curvature  (see  {B.69))  can  be  written  as 


Up 
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J^k*'^E{k*)dk*  2 

'jQ^E(k*)dk*  Jo 


-t-  6j)  cos^ 6d9. 


(2.17) 


As  seen  in  (2.16)  and  (2.17),  the  statistics  of  shock  front  displacement  and  its 
curvature  are  dependent  on  the  shape  of  the  upstream  spectrum,  E{k*).  The 
shock  wave  displacement  has  significant  contributions  from  large  scale  turbulence. 
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while  its  curvature  is  scaled  with  the  inverse  of  upstream  turbulence  microscale 

J  k^E{k)dk/  J  E{k)dk  ~  1/A2). 

Statistics  of  shock  front  distortion  are  obtained  by  numerically  integrating  (5.65), 
(2.16)  and  (2.17)  using  the  energy  spectrum  in  (2.1).  Figure  2.12  presents  the  rms 
Vcdues  of 


Note  that  the  statistics  of  the  shock  front  distortions  are  seeded  with  the  upstream 
turbulence  intensity.  As  the  mean  upstream  Mach  number  increases,  the  scaled 
rms  values  are  found  to  decrease. 

Considering  the  time  dependence  of  the  upstream  velocity  fluctuations  at  a  fixed 
point,  LI  A  can  predict  the  local  fluctuating  shock  front  speed,  (f.  Expression 
(5.27)  can  be  rewritten  as 

dui  =  dui  exp  [i  (fci(xi  —  U\t)  +  Ar2X2  +  ^3^3)]  •  (2.18) 

Since  the  velocity  fluctuations  at  the  shock  front  (xj  =  0)  vary  not  only  in  the 
transverse  direction  but  also  in  time,  the  argument  of  the  exponential  function 
in  (2.18)  becomes  i  (  — kiFit  +  k-2X2  4- 1:3x3)  with  xj  =  0.  Therefore,  expression 
(5.66)  for  the  local  shock  displacement  can  be  rewritten  as 

i  =  T^^exp[i(k;,  -x/i  -  k^Uit  +  ^sj,)] ,  (2.19) 

ikh 


where  is  the  magnitude  of  kf^  =  (A:2i^'3)i  =  (^21^:3).  From  (2.19), 

one  obtains  the  relation  between  the  local  fluctuating  shock  front  speed,  and 
the  local  shock  front  inclination  angle  in  the  transverse  direction  (x2  direction), 
<72  =  ^,2i  as 


kiUi 

^2 


^2- 


(2.20) 


Using  (5.57)  and  (5.58),  (2.20)  can  be  expressed  as 


tt 


•V 


f<n 


hlcx]y[i{kf,  ■  X/,  +  -^sh)!- 


(2.21) 
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Figure  2.13  compares  the  local  upstream  velocity  fluctuation  and  the  local  shock 
front  speed  for  upstream  waves  incident  on  a  shock  wave  of  =  1.2  at  different 
angles.  The  magnitude  of  the  shock  speed  is  comparable  to  the  upstream  fluctuation 
velocity.  For  waves  whose  incident  angles  are  smaller  than  the  critical  angle  of 
incidence,  =  36.4°,  the  local  shock  front  speed  lags  by  a  phase  angle  between  0 
to  90°  with  respect  to  the  upstream  velocity  fluctuation;  for  waves  whose  incident 
angles  are  larger  than  the  critical  angle,  the  shock  front  speed  is  in  phase  with 
the  velocity  fluctuation.  The  local  shock  front  speed  which  is  in  phase  with  the 
upstream  fluctuation  velocity  attenuates  the  fluctuations  in  the  effective  upstream 
Mach  number,  =  {U\  +  uj  — 

Figures  2.14(a)  and  (b)  show  the  dependence  of  the  rms  fluctuating  shock  front 
speed  on  the  mean  upstream  Mach  number.  The  shock  front  speed  exceeds  the 
upstream  fluctuation  velocity  for  <  1.25.  For  weak  shock  waves,  the  shock 
front  speed  is  close  to  the  upstream  fluctuation  velocity  resulting  in  approximately 
uniform  effective  upstream  Mach  number. 
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Figure  2.1.  the  coordinate  system  used  in  the  linear  analyses. 


Figure  2.2.  Predictions  of  transverse  vorticity  amplification:  (a)  1  <  <  10, 

(fc)  1  <  m['  <  2.  - T.IA, - RDT. 
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Figure  2.3.  Predictions  of  solenoidal  velocity  component  amplification:  (a)  1  < 
m['  <  10,  (h)  1  <  <  2.  -  «i(LIA). - - ui(RDT), 
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Deray  of  velocity  coniponenls  behind  the  shock  wave  for  expor.  ntial  fall-off  upstream  spectrum  (2.1), 
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Figure  2.7.  Decay  of  stream  wise  velocity  fluctuation  behind  the  shock  wave  for  a  von  Karman  upstream  spectrum 
(2.2),  with  mF  =  1.2:  - - U2^u^. 


cf 


CURE  2.8.  Frequency  spectrum  nmplificnlion  ratio  for  a  von  Karman  up 
•earn  spectrum,  (2.2);  (a)  RDT,  (t)  LIA. -  =  11. 


Figure  2.10.  Predictions  of  changes  in  Taylor  microscales:  (a)  1  <  Mf  <  10, (6) 

1  <  <  2.  - Ai(LIA), - A2&A3(LIA), - Ai(RDT),  A2&A3(RD1 
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Figure  2.11.  Acoustic  pressure  fluctuation  downstream  of  the  shock  wave:  (a) 
normalized  with  the  downstream  mean  pressure,  (6)  normalized  with  the  upstream 
mean  pressure.  - near  field, - far  field. 


Figure  2.13.  Comparison  of  local  upstream  velocity  fluctuation  and  the  local  shock  front  speed  for  various  ang 
of  incidence  (ref.  (2.18)):  - ui, - ^,<(«  =  15°), - =  30°), .  ^,<(0  =  45°), - ^,<(<?  =  60°) 


CHAPTER  3 


NUMERICAL  METHOD 

The  time-dependent  Navier-Stokes  equations  for  a  compressible  fluid  were  solved 
directly.  All  the  relevant  turbulence  scales  are  resolved  without  a  turbulence  model, 
and  the  shock  wave  structure  is  resolved  as  a  solution  of  the  Navier-Stokes  equa¬ 
tions  without  introducing  the  techniques  of  shock-fitting  or  shock-capturing.  The 
shock  structure  is  adequately  represented  by  the  Navier-Stokes  equations  for  Mach 
numbers  less  than  2.0  [Sherman  1955].  Except  for  monatomic  gases,  how¬ 
ever,  the  thickness  of  the  shock  wave  as  a  solution  of  the  Navier-Stokes  equations 
is  underpredicted  even  for  <  2.0,  because  the  rotational  energy  mode  is  not 
in  equilibrium  inside  the  shock  wave  [Lumpkin  1990].  This  chapter  describes  the 
governing  equations  and  the  numerictd  method  used  to  simulate  shock/turbulence 
interaction,  where  the  flow  is  assumed  to  be  periodic  in  transverse  directions  {x2 
and  13).  This  chapter  also  includes  validations  of  our  numerical  schemes.  The  code 
is  written  in  the  VECTORAL  language  [Wray  1988]  and  implemented  on  a  Cray 
Y-MP/832  at  NASA-Ames  Resea-ch  Center. 


3.1  Governing  Equations 

The  conservation  laws  for  mass,  momentum,  and  energy  are  [Anderson,  Tannehill 
and  Fletcher  1984]; 


dp*  d{p*u*)  _ 
dt*  dx* 

d{p*u*)  9{p*u*u*  ^  p*6tj)  dr*- 

dt*  dx*  ~  ^ 

dE*^  d[{E*^  +  p*)u*]  dQ*  diu*r*^) 

dt*  ^  dx*  dx*  ^  dx* 


(3.1) 

(3.2) 


(3.3) 


where  superscript  *  indicates  a  dimensional  quantity,  p*  is  the  density,  u*  the 
velocity  components,  p*  the  pressure,  t*^  the  shear  stress  tensor,  and  Q*  is  the  heat 
flux  vector.  Ej  is  the  total  energy,  defined  by: 
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(3.4) 


u*u* 

where  e*  is  the  internal  energy  per  unit  mass. 

We  assume  the  fluid  to  be  a  perfect  gas  satisfying 

p*  =  p*R*T\  (3.5) 

where  R*  is  the  gas  constant  and  T*  the  temperature.  We  assume  a  Newtonian 
fluid  and  use  Stokes  hypothesis  and  Fourier  law  of  heat  conduction,  so  that  the 
constitutive  equations  for  r-*  and  Q*  are: 


M  ( 


du)  2dul 


dx*-  dx*  3  3zi 

J  *  * 


Q*  =  -k 


^dT* 

ax- 


where  p*  is  the  molecular  viscosity  and  k*  the  thermal  conductivity. 
The  flow  variables  are  non-dimensionalized  as  follows 


u 


Ui  = 


P 


P  = 


Po<:f 


T  = 


T* 


(7-i)r* 


p 

p  =  -i 
Po 


t  = 


t*i 


Xi  = 


LI' 


(3.6) 

(3.7) 


(3.8a) 


(3.86) 


where  subscript  o  represents  the  mean  upstream  value,  c*  is  the  speed  of  sound, 
and  Lq  is  an  arbitrary  reference  length  scale.  The  value  of  the  specific  heat  ratio  7 
defined  as 

7  =  Cp/c*  (3.9) 

is  taken  as  1.40  in  this  work. 

The  non-dimensiontd  equations  for  continuity,  momentum,  and  energy  are 


^  ^  d{pu^) 
dt  dxi 


(3.10) 
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(3.11) 


djpuj)  _  djpu^uj  +  pS^j)  dTjj 
dt  dxj  dxj 

dEp  _  d{{Ex  +  p)ut]  dQi  dujTjj 

dt  dxi  dxf  ^  dxj 

with  constitutive  relations 


^  _  2dui, 

Re  dxj  dxi  3  dxf.  ’ 


Q^  =  - 


p  dr 

PrRe  dx^ 


where  Re  and  Pr  are  defined  as 


Re  = 


PlclLl 

P*o 


Pr  = 


We  assume  the  Prandtl  number  to  be  constant  equal  to  0.70  and 
cosity  to  follow  the  power  law: 


Po 


or  p  =  ((7-l)r]", 


where  n  =  0.76. 

If  we  assume  constant  specific  heats  and  set  e*  =  0  at  = 
e*  =  c^T* .  The  non-dimensional  form  of  the  perfect  gas  law  is 


P  = 


(7-1)  ^ 

- pT. 

7 


(3.12) 

(3.13) 

(3.14) 

(3.15a) 

(3.156) 
assume  the  vis- 

(3.16) 

0,  we  can  write 

(3.17) 
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3.2  Numerical  Schemes 


This  section  describes  the  numerical  schemes  used:  time  advancement  scheme, 
approximation  of  spatial  derivatives,  special  treatment  of  convective  terms  for  nu¬ 
merical  stability,  method  of  generating  inflow  turbulence,  and  initial  and  boundary 
conditions. 

3.2.1  Time  Advancement 

An  explicit  time-advancement  method  is  used.  The  variables  {p,pui,  Ej')  are  ad¬ 
vanced  using  a  three-step  compact-storage  third-order  Runge-Kutta  scheme  [Wray 
1986).  This  scheme,  when  applied  to  dy/dt  =  has  the  following  three  sub¬ 

steps: 


=  V"  +  i/(!,",(»)A(  +  j/(!,“,(“)A(, 

where  -f  and  t'*  =  <”  -f  jAt. 

The  time  step  is  computed  from  the  following  formula: 

lAlmax’ 


(3.18a) 

(3.186) 

(3.18c) 


(3.19) 


where  C FL  is  the  Courant-Friedrichs-Lewy  number.  The  subscript  max  refers  to  the 
maximum  over  all  grid  points.  The  maximum  CFL  number  for  stability  is  fixed  by 
the  time-advance  method.  For  linear  equations  the  limit  is  \/^  for  the  third-order 
Runge-Kutta  method  described  above.  For  the  three-dimensional  Navier-Stokes 
equations,  we  have  used: 


(W\\  +  c  M  |m3I  +  c\ 

\  Axi  Aar2  Ax3  / 


(3.20) 
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where  |ui|  is  the  absolute  velocity  in  the  direction  and  Ax^  is  the  grid  spacing. 
3.2.2  Evaluation  of  Spatial  Derivatives 


A  family  of  high-order  modified  Fade  schemes  has  been  derived  by  Lele  [1990] 
with  spectral-like  resolution  characteristics  as  well  as  high-order  formal  accuracy. 
In  this  work,  we  used  such  a  scheme  for  both  the  first  and  second  derivatives. 

The  first  derivative  is  given  by: 


/  .  /  .  /  .  y>+i  -  J/y-i  .  yj+2  -  yj-2 

yj-l  +  <■!»>  +  Vj+l  =  02  — 2Ai —  +  4AJ — ■ 


(3.21a) 


We  can  obtain  y'j  by  solving  a  tridiagoncil  system  of  equations.  A  family  of  fourth 
order  schemes  is  obtained  if  we  choose 


“2 


2  4"  Adi 

3 


03  = 


A  —  a\ 
3 


(3.216) 


For  aj  =  4,  the  conventional  Fade  scheme  is  recovered,  while  with  C]  =  3  and 
(3.216)  we  have  a  sixth-order  scheme.  Similarly,  we  can  write  for  the  second  deriva¬ 
tive: 


n  I  n  ,  n 
yj-l  +aiyj  +  yj^i  =  02 


yj+i  -  2yj  +  y>-i  Vj+2  -  2yj  +  yj-2 


Ai2 


+  03 


4A5 


(3.22a) 


4ai  —  4  10  —  oj 

<■2  =  -  3—  a,  = 


(3.226) 


The  choice  of  aj  =  10  recovers  the  usual  Fade  scheme,  and  for  aj  =  11/2  with 
(3,226)  it  is  sixth-order  accurate. 

In  this  work,  we  used  the  sixth-order  schemes  both  for  the  first  and  second 
derivatives.  The  resolution  characteristics  of  the  conventional  Fade  scheme  are 
improved  with  negligible  increase  of  operations. 


The  diffusion  terms  in  the  governing  equations  require  evaluation  of  successive 
derivatives,  for  example 


When  a  spectral  method  is  used,  there  is  no  loss  of  accuracy  when  these  are  com¬ 
puted  by  two  applications  of  a  first-derivative  operator.  With  finite  difference 
methods,  however,  two  applications  of  a  first  derivative  results  in  a  less  accurate 
representation  of  the  derivative  at  high  wave  numbers  as  compared  to  one  applica¬ 
tion  of  a  second-derivative  operator,  for  the  modified  wave  number  goes  to  zero  for 
the  first  derivative  at  high  wave  numbers  (Lele  1990].  To  eliminate  this  inaccuracy, 
we  expanded  all  diffusion  terms  into  two  terms  (non-conservative  formulation): 


d^ui  dfJ.  dui 
^  ^  dxi ' 


(3.24) 


and  use  the  formulae  (3.21)  and  (3.22)  to  approximate  the  first  and  the  second 
derivatives,  respectively. 

To  simulate  shock/turbulence  interaction,  a  shock  wave  is  placed  near  the  center 
of  the  computational  domain  (Figure  3.1)  and  numerical  simulations  are  performed 
in  a  frame  moving  with  the  shock  wave.  Since  the  resolution  requirements  for  a 
shock  wave  are  far  more  restrictive  than  those  for  turbulence,  a  non-uniform  grid 
is  used  in  the  streamwise  (or  xj)  direction  to  resolve  the  shock  wave  structure. 
The  following  mapping  from  the  uniform  computational  grid  to  the  non-uniform 
physical  grid  is  used  to  concentrate  points  near  the  region  occupied  by  the  shock 
wave: 


Li 


s  —  ^^erf(cis) 

1  -  v^'^erf(ci/2) 


(3.25) 


where  Li  is  the  length  of  the  computational  box  in  the  streamwise  direction,  s  is  the 
coordinate  in  the  computational  grid  ranging  between  (  — 2,2),  and  61  and  cj  are 
the  stretching  parameters  controlling  grid  stretching  ratio  and  the  grid  stretching 
rate  respectively.  The  grid  stretching  ratio,  (Axi)max/(^a^l)min  —  1/(1  “  ^l)i  was 
chosen  to  be  between  5  and  20.  Higher  grid  stretching  is  required  for  a  larger 
separation  between  the  shock  wave  thickness  and  turbulence  length  scales.  W^ith 
larger  cj,  the  rate  of  grid  stretching  is  faster.  In  this  work  cj  is  chosen  to  be  5.  If 
we  define  the  metric  quantities: 
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3.2.3  Nonlinear  Numerical  Stability 


Compressible  flow  simulations  with  conservative  formulation  are  especially  prone 
to  aliasing  errors  because  evaluations  of  velocity  and  temperature  from  the  conser¬ 
vative  variables  involve  the  division  operation  which  has  no  clear  interpretation  in 
the  Fourier  space.  A  possible  way  of  conducting  alias-free  simulation  of  compress¬ 
ible  turbulence  is  to  solve  for  the  specific  volume  in  a  mass  conservation  equation. 
We  discuss  this  alternative  in  Appendix  C. 

Feiereisen  et  al.  [1981]  have  noted  that  the  use  of  a  special  form  of  the  convection 
term  with  a  symmetric  differencing  in  space  ensures  conservation  of  total  energy 
in  the  inviscid  limit.  Blaisdell  et  al.  [1990]  have  explicitly  shown  that  this  special 
form  helps  to  control  aliasing  errors.  For  the  numerical  stability,  we  evaluate  the 
convection  term  in  the  momentum  equations  in  the  following  special  form: 


dpuiUj  1  Idpu^Uj  duf  dpuj' 

dxj  2  dxj  ^  ^  ^  Qxj  ^  ‘  dxj 


(3.30) 
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In  a  spatially  evolving  simulation,  spurious  numerical  waves  are  generated  at  the 
inflow  as  soon  as  disturbances  encounter  the  outflow  boundary  [Buell  and  Huerre, 
1988].  These  numerical  waves  have  a  wave  length  of  twice  the  grid  size  and  continue 
to  grow  in  time.  To  remove  these  spurious  waves,  localized  filtering  near  the  inflow 
plane  is  performed  using  the  following  scheme: 


Vi-l  +  aiVi  +  =  0-2  (y,_i  +  2yi  +  +  03  (y,_2  -  2yi  +  yi+2),  (3.31) 

where  y  and  y^  are  the  unfiltered  and  the  filtered  quantities  respectively.  We  choose 


02 


ai  +  2 


03 


2  —  0} 
16  ■ 


The  unfiltered  quantity  is  recovered  for  aj  =  2.  This  filtering  operation  exactly  re¬ 
moves  waves  of  twice  the  grid  size.  The  filter  transfer  function  in  the  Fourier  space, 
defined  as  the  ratio  of  the  Fourier  transforms  of  y^  and  y,  is  shown  in  Figure  3.2. 
For  aj  close  to  2,  filtering  operation  is  more  localized  at  high  wave  numbers.  Since 
the  spurious  numerical  waves  travel  mainly  in  the  streamwise  direction,  appearing 
first  at  the  inflow  boundary  with  wave  length  of  twice  the  grid  size,  the  filtering 
operation  is  performed  only  in  the  streamwise  direction  near  the  inflow  boundary 
over  less  than  a  tenth  of  the  computationcd  domain  with  aj  =  2.01. 


3.2.4  Inflow/Outflow  Boundary  Conditions 


Many  of  the  existing  boundary  conditions  for  the  compressible  Navier-Stokes 
equations  are  based  on  the  concept  of  the  characteristics  along  which  information 
travels.  The  number  of  boundary  conditions  required  at  a  point  on  a  boundary 
varies  with  the  flow  conditions  at  that  point:  in  general,  the  number  of  boundary 
conditions  which  must  be  specified  at  a  point  on  a  boundary  is  equal  to  the  number 
of  incoming  waves  (from  outside  of  the  domeiin)  at  that  point  [Thompson  1987]. 

In  this  work,  the  inflow  is  kept  supersonic  so  that  we  can  specify  all  the  flow 
variables.  This  is  because  supersonic  inflow  guarantees  all  the  information  to  be 
incoming.  Mean  values  of  velocity,  pressure,  and  density  are  set  to  be  constant  over 
the  inflow  plane.  The  turbulence  velocity  signal  generated  at  the  inflow  boundary 
is  designed  to  be  isotropic  with  a  prescribed  spectrum  with  no  fluctuations  in 
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pressure  and  density.  Appendix  D  gives  a  detailed  description  and  validations  of 
the  procedure  for  generation  of  inflow  turbulen  :e. 

To  generate  inflow  turbulence,  we  use  the  following  three-dimensional  energy 
spectrum  function: 


E{k)  =  16^^(A)4exp[-2(A)2]^  (3.32) 

where  Uq  is  the  rms  turbulence  intensity  and  ko  is  the  most  energetic  wave  number. 
This  spectrum  has  the  following  properties: 


_2  foo  o 

^  =  E{k)dk  =  -ul 

(3.33) 

€  =  2fj,  f  k^E{k)dk  =  —fiu^k^ 

Jo  4 

(3.34) 

u  uko 

(3.35) 

where  u  is  the  kinematic  viscosity,  and  A  is  the  longitudinal  Taylor  microscede. 

Outflow  is  subsonic  in  the  mean  sense,  which  requires  special  attention  in  ap¬ 
plying  boundary  conditions,  since  one  of  the  characteristics  is  incoming  (from  out¬ 
side).  We  can  determine  what  information  the  outgoing  characteristics  carry  from 
the  solution  inside  the  computational  domain.  But  difficulty  arises  in  determin¬ 
ing  contributions  from  the  incoming  characteristics.  Since  we  are  limited  to  the 
information  inside  the  computational  domain,  the  information  that  can  be  used  to 
determine  the  effect  of  incoming  characteristics  is  incomplete,  and  we  need  to  make 
some  assumptions. 

Thompson’s  boundary  condition  [1987]  is  derived  from  the  assumption  that  the 
unknown  incoming  data  have  no  effect  on  the  flow  variables  at  the  boundaries. 
Poinsot  et  al.  [1990]  extended  Thompson’s  boundary  condition  to  include  viscous 
effects  at  the  boundaries.  Giles’  [1990]  boundary  condition  is  designed  in  the  Fourier 
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space  so  that  small  amplitude  waves  may  leave  the  computational  domain  without 
reflection. 

In  this  work,  we  tested  the  boundary  conditions  by  Thompson  [1987],  Poinsot 
et  al.  [1990],  and  Giles  [1990].  We  found  that  they  were  comparable  in  suppressing 
numerical  reflections  at  the  boundary,  with  Giles’  boundary  condition  being  supe¬ 
rior  to  the  other  two.  However,  this  improvement  had  virtually  no  effect  on  the 
turbulence  statistics  downstream  of  the  shock  wave  (for  details,  see  Appendix  G). 
For  simplicity  of  implementation,  we  used  Thompson’s  method  in  most  of  this  work. 
In  the  remaining  part  of  this  subsection,  we  show  the  derivation  and  application  of 
Thompson’s  boundary  condition. 

The  basic  idea  is  to  consider  the  characteristic  form  of  the  Euler  equations  at 
the  outflow  boundary.  Outgoing  characteristics  use  information  from  within  the 
computational  domain,  and  can  be  computed  with  no  difficulty.  Incoming  char¬ 
acteristics  are  handled  by  setting  the  time  derivative  of  their  amplitude  equal  to 
zero,  thus  ensuring  that  no  waves  enter  the  domain  during  the  simulation,  giving 
the  boundary  conditions  a  non-reflecting  character. 

VV’e  begin  by  writing  the  Euler  equations  in  terms  of  the  conservative  variables 


dt  ^  dx  I 


=  (RHS), 


(3.36) 


where 


F  (^pui,pu]  -t-  p,puiU2,puiU2,{ET  +  P)“l) 


(3.37) 


We  are  concerned  here  only  with  the  derivatives  in  the  ij  direction.  Derivatives 
in  T'}  and  ^3,  and  visccjus  terms,  are  evaluated  at  the  boundary  using  information 
from  ^he  previous  substep,  and  are  included  in  the  right  hand  side  (RHS).  The 
flux  .lacobian  of  F  is  more  easily  derived  if  we  work  with  the  non-conservative  flow 
variables  U  -  (p,  u j ,  U2,  U3,  p)^.  Setting  A  =  dFjdQ  (t.e.  =  dFj/dQj)  and 

R  -  OQ'dJF  we  have 


dV  ,  dU 

—  •F  .4  — 

()i  Ox  J 


(RHS) 


(3.38) 
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and 
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(3.39) 


Now  A  can  be  diagonalized,  A  =  T  ^AT,  where  the  elements  of  the  diagonal 
matrix  A  are  (u^  —  c,iii  +  c,ui,ui,ui).  Equation  (3.39)  can  now  be  written  as 


(3.40) 


This  is  the  relation  that  is  imposed  at  the  boundary  to  calculate  dF/dxi  in  (3.36). 

The  quantity  in  the  parenthesis  in  (3.40)  is  a  vector.  The  sign  of  each  eigenvalue 
in  A  is  used  to  determine  the  course  of  action  for  each  element  in  the  vector.  If 
the  characteristic  velocity  is  directed  out  of  the  computational  domain  (positive 
eigenvalue  at  the  outflow  boundary),  then  the  quantity  is  calculated  as  it  stands 
using  a  one-sided  difference.  On  the  other  hand,  if  the  characteristic  is  directed 
inwards  then  the  element  is  set  to  zero.  This  gives  the  non-reflecting  character  of 
the  boundary  condition  for  waves  at  close  to  normzd  incidence.  The  final  step  is  to 
premultiply  the  vector  by  the  matrices  T~^  and  R.  The  various  matrices  required 
in  the  computation  are: 
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3.2.5  Initial  Conditions 

To  simulate  shock  turbulence  interaction,  we  initialize  the  field  by  superposing 
isotropic  homogeneous  turbulence  on  the  corresponding  stationary  laminar  shock 
wave  profile.  This  presumably  reduces  the  time  required  for  the  shock  wave  and 
turbulence  to  reach  statistically  stationary  state.  The  turbulence  field  used  in 
the  initial  condition  has  the  same  energy  spectrum  as  turbulence  generated  at 
the  inflow  boundary  during  the  computation.  Rogallo’s  scheme  [1981]  is  used  for 
turbulence  initialization,  and  attention  is  given  to  ensure  smoothness  of  turbulence 
signal  at  inflow  boundary.  We  solved  the  Navier-Stokes  equations  to  get  a  laminar 
^  shock  wave,  and  used  this  solution  to  initialize  a  planar  laminar  shock  wave  in  the 

computational  domain. 

3.2.0  Statistical  Averages 

Two  different  statistical  averages  are  used;  the  conventional  ensemble  average 
and  the  Favre  (or  density-weighted)  ensemble  average.  The  ensemble  average  of 
a  quantity  is  defined  as  the  average  over  time  and  over  the 

(homogeneous)  transverse  directions.  The  en.semble-averaged  equations  in  com¬ 
pressible  turbulence  are  usually  too  complex.  The  averaged  equations  are  simpler 
using  density-weighted  averages  suggested  by  Favre  [1965a,  1965b].  The  density- 
weighted  average  y{x])  is  defined  as 

.V(^l)--^  (3.44) 

r 

The  fluctuations  from  the  ensemlile  aver.uge  and  from  the  Favre  average  are  defined 
as 

y'  y  .v(>i).  y”  ^  y  -  y{x\). 
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The  ensemble  average  of  y"  is  not  zero;  instead, 


y 
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=  y- 


py  +  p'y’ 
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p'y' 

p 


(3.45) 


3.3  Validations 

In  this  section,  we  provide  a  validation  of  the  computer  code  by  separate  com¬ 
putations  of  some  of  the  components  of  the  shock/turbulence  interaction  problem. 
For  a  direct  numerical  simulation  of  spatially  evolving  turbulence  interacting  with  a 
shock  wave,  the  shock  wave  needs  to  be  well  resolved,  spati2dly  evolving  turbulence 
must  be  properly  simulated,  and  the  interaction  of  turbulence  with  a  shock  wave 
must  be  accurately  predicted.  For  this  purpose  there  are  at  least  three  independent 
categories  to  be  validated:  (1)  simulation  of  spatially  evolving  turbulence  (without 
a  shock  wave),  (2)  resolution  of  shock  wave  structure,  and  (3)  interaction  of  flow 
inhomogeneity  with  a  shock  wave. 

The  validation  of  the  first  category  is  given  in  Appendix  D  by  comparing  the 
simulation  results  of  spatially  evolving  turbulence  with  the  experimental  data. 

Validations  of  the  other  two  categories  are  performed  in  the  present  section. 
The  capability  of  resolving  the  shock  wave  structure  is  validated  by  solving  the 
one-dimensional  Navier-Stokes  equations,  investigating  the  shock  wave  profile  qual¬ 
itatively  and  comparing  the  shock  wave  thickness  with  a  theoretical  prediction. 
Subsequently,  interaction  of  a  flow  inhoniogeneity  with  a  shock  wave  is  checked  by 
investigating  vorticity  production  through  the  baroclinic  torque  and  comparing  the 
resulting  circulation  with  a  theoretical  prediction. 
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3.3.1  Laminar  Shock  Wave 


The  shock  wave  structure  was  resolved  by  solving  the  Navier-Stokes  equations. 
Prediction  of  profiles  across  the  shock  wave  was  found  rehable  up  to  an  upstream 
Mach  number  of  2.0  [Sherman  1955].  Beyond  this  Mach  number,  the  thermal  equi¬ 
librium  assumption  is  no  longer  valid,  preventing  use  of  the  Navier-Stokes  equations 
to  resolve  the  shock  wave  structure.  Except  for  monatomic  gases,  however,  the 
thickness  of  the  shock  wave  as  a  solution  of  the  Navier-Stokes  equations  is  found 
to  be  underpredicted  even  for  <  2.0  [Lumpkin  1990]. 

In  the  following,  we  test  the  ability  of  our  numerical  scheme  to  resolve  a  normcd 
shock  wave  in  steady  flow  and  estimate  the  number  of  grid  points  required  across 
a  shock  wave  to  properly  resolve  its  structure.  We  computed  the  normal  shock 
wave  for  the  upstream  Mach  number,  —  1.2,  where  superscript  ^  denotes  an 
upstream  value.  Calculations  are  performed  in  a  frame  fixed  on  the  shock  wave 
where  the  inflow  is  supersonic  and  the  outflow'  subsonic  with  uniformly  distributed 
201  grid  points. 

The  initial  conditions  of  density,  pressure,  and  velocity  are  given  by  the  following 
expression,  which  satisfies  the  Rankine-Hugoniot  relations  across  the  shock  wave; 

y{xi)  =  +  (y^  -  y^)  tanh(^i-^^),  (3.36) 

where  superscript  ^  represents  a  downstream  value,  y  is  one  of  the  flow  variables, 
and  Xj  and  6*  represent  the  shock  center  position  and  the  initial  shock  thickness 
parameter,  respectively.  The  flow  variables  are  fixed  at  the  supersonic  inflow  and 
Thompson’s  non-reflecting  boundary  conditions  are  applied  at  the  subsonic  outflow 
boundary.  As  the  Navier-Stokes  equations  are  advanced,  initial  profiles  relax  or 
steepen  into  equilibrium  profiles. 

The  profiles  of  velocity,  pressure,  temperature,  dilatation,  and  entropy  across  the 
shock  wave  are  shown  in  Figures  3.3(a-e).  All  the  Rankine-Hugoniot  conditions  are 
satisfied  without  any  spurious  oscillations  across  the  shock  wave.  The  dilatation 
profile  shows  that  the  shock  wave  is  well  resolved  without  any  numerical  problems. 
The  entropy  profile  shows  a  local  peak  inside  the  compression  zone  as  well  as  a  net 
increase  across  the  zone.  As  flow  passes  through  a  shock  wave,  flow  kinetic  energy 
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is  converted  into  internal  energy  by  viscous  dissipation,  which  leads  to  entropy 
production;  the  resulting  temperature  gradient  leads  to  further  irreversible  entropy 
production  and  reversible  entropy  flux.  Irreversibility  from  viscous  dissipation  and 
temperature  gradient  are  responsible  for  the  overall  increase  of  entropy,  while  the 
reversible  entropy  flux  leads  to  the  peak  inside  the  shock  [Lagerstrom  1964].  Figure 
3.3(f)  shows  the  budget  of  terms  in  the  entropy  transport  equation, 

Jb  /dry 

^  ^  dx\  dxi  \T  dxi )  ^  ZT  \dxi  )  ^  T^\dx\) 

' - V - '  ' - S/ - '  ' - V - '  ' - V.. - ' 

(I)  (H)  (III)  (IV) 

[Thompson  1984],  where  the  convection  (I)  is  balanced  by  the  sum  of  the  reversible 
entropy  flux  (II)  and  the  entropy  production  by  viscous  dissipation  (III)  and  irre¬ 
versible  heat  transfer  (IV).  The  figure  verifies  that  the  entropy  decrease  inside  the 
shock  wave  is  due  to  the  reversible  entropy  flux  (II).  The  primary  source  of  the  net 
entropy  increase  across  the  shock  wave  is  found  to  be  the  viscous  dissipation  of  the 
flow  kinetic  energy  (III). 

The  shock  wave  thickness,  defined  as 


it/f  -  t/f  I 

Wild. 


Imax 


(3.37) 


is  estimated  to  be 


7  +  1  \3  Pr  J  Re{M*  -  1) 


(3.38) 


for  very  weak  shocks  [Shapiro  1953],  where  Re  is  defined  as  (3.15a)  and 


1  + 
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The  shock  wave  thickness  in  the  simulation  is  within  7  percent  of  the  estimate  in 
(3.38)  when  the  grid  spacing  in  the  shock  zone  is  less  than  a  third  of  the  shock 
wave  thickness,  that  is,  Aij  <  36,. 

3.3.2  Thermal  Inhomogeneity  Interacting  with  a  Shock  Wave 

A  second  validation  study  was  performed  to  check  the  accuracy  of  the  scheme 
in  the  prediction  of  the  interaction  of  time-dependent  disturbances  with  a  shock 
wave.  This  study  also  provided  an  evaluation  of  Thompson’s  boundary  conditions 
at  the  unsteady  subsonic  outflow  boundary. 

The  interaction  of  a  thermal  inhomogeneity  with  a  shock  wave  with  =  1.20 
was  simulated.  Inhomogeneity  at  the  inflow  boundary  had  a  circular  shape  in  the 
(t-X2)  plane  as 


T(x2,i)  = 

7  -  1 


+  aj  exp 


(3.39) 


where  aj  is  the  relative  amplitude  of  the  disturbance,  f  j  the  time  at  which  the  cen¬ 
ter  of  the  disturbance  passes  through  the  inflow  plane,  bj  the  length  scale  of  the 
disturbance,  and  nj'  the  parameter  controlling  the  sharpness  of  the  disturbance’s 
edge.  The  simulation  used  201  x  101  uniform  meshes.  The  values  used  in  the  simu¬ 
lation  are  aj  =  1/10,67’  —  £'2/4,  and  nj  =  10.  {L2  is  the  size  of  the  computational 
box  in  the  X2  direction.)  Pressure  at  the  inflow  plane  is  kept  uniform  and  constant. 
Thompson’s  non-reflecting  boundary  condition  is  used  at  the  outflow  boundary  and 
periodic  boundary  conditions  are  used  in  the  transverse  direction. 

Figures  3.4(a-c)  show  a  time  sequence  of  the  temperature  field.  Upstream  of  the 
shock  wave  the  shape  of  the  iiihomogeneity  is  circular,  and  after  the  interaction 
the  shape  is  changed  into  an  ellipse.  The  predicted  aspect  ratio  of  this  ellipse 
is  within  one  percent  of  the  exact  value,  equal  to  the  density  ratio  across  the 
shock  wave,  j .  In  Figures  3..')(a-c),  a  t’me  sequence  of  the  vorticity  field  is 
shown.  Through  the  interaction  of  the  thermal  inhomogeneity  with  the  shock  wave, 
vorticity  is  produced  by  the  baroclinic  torque.  Figure  3.5  clearly  shows  counter¬ 
rotating  vortices  generated  during  the  interaction. 
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Picone  et  al.  [1985]  have  derived  an  expression  for  the  circulation  of  the  vor¬ 
tices  generated  by  the  baroclinic  torque,  assuming  the  shock  wave  to  be  a  planar 
discontinuity  and  the  thermal  inhomogeneity  to  be  a  circular  discontinuity.  Their 
expression  for  the  circulation  in  the  upper  half  plane  is 

r  ~  (<  -  t'f)  Ml  +  “t).  (3.40) 

In  Figure  3.6,  the  history  of  circulation  in  the  upper  half  plane  is  shown.  We  see 
that  circulation  peaks  near  the  end  of  the  interaction  and  decays  through  viscous 
diffusion  as  vortices  flow  downstream.  This  peak  strength  compares  favorably  with 
the  estimation  of  (3.40)  to  within  5%,  thus  confirming  the  ability  of  the  scheme 
to  predict  the  shock/disturbance  interaction  correctly.  More  rapid  decrease  in 
circulation  after  t  =  20  is  due  to  the  primary  vortices  leaving  the  domain.  There  are 
treuling  vortices  with  circulations  of  opposite  signs  to  the  main  vortices  generated 
through  the  interaction  (see  Figure  3.7).  The  trailing  vortices  are  prodi-.c  d  by  the 
relaxation  of  the  curvature  in  the  shock  wave,  and  their  strengths  are  an  order  of 
magnitude  weaker  than  the  main  vortices. 

To  check  the  accuracy  of  Thompson’s  non-reflecting  boundary  condition  in  un¬ 
steady  flows,  we  examined  the  pressure  field  at  times  of  entrance,  interaction,  and 
exit  of  thermal  inhomogeneity  in  Figures  3.8(a-c).  The  disturbance  pressure  at  the 
outflow  boundary  is  not  noticed  in  these  plots.  The  level  of  disturbance  pressure  af¬ 
ter  the  exit  of  the  inhomogeneity  from  the  computational  domain  is  about  0(10”^) 
compared  to  the  mean  pressure. 
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3.2.  The  transfer  function  of  filtering  in  t! 
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Figure  3.3.  Continued.  For  Figure  3.3(/):  -  convection,  entropy 

flux, - production  by  viscous  dissipation,  — ■ —  production  by  irreversible  heat 

transfer. 
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Figure  3.8.  Pressure  fields  from  passage  of  a  thermal  inhomogeneity  through  £ 
shock  (The  contour  increments  are  7Ap  =  0.05  starting  from  the  mean  upstrean 
pressure  7p  =  1.0).  (a)  t  =  3.00,  (b)  t  =  11.0,  (c)  <  =  22.2. 


CHAPTER  4 


DIRECT  NUMERICAL  SIMULATION 


In  this  chapter  we  discuss  the  results  from  direct  numerical  simulations  of  the 
interaction  between  isotropic  turbulence  and  a  normal  shock  wave. 

The  simulations  are  conducted  in  a  reference  frame  fixed  with  respect  to  the 
mean  shock  position  so  that  long-time  statistical  averages  of  turbulence  quantities 
could  be  obtained.  In  this  frame  of  reference,  the  mean  flow  approaches  the  shock 
wave  with  a  supersonic  speed  and  leaves  with  a  subsonic  speed.  Inflow  turbulence  is 
generated  using  the  method  described  in  Appendix  D,  and  the  pressure  and  density 
are  kept  constant  and  uniform  in  the  inflow  plane.  Fluctuations  in  pressure  and 
density  naturally  evolve  as  the  flow  approaches  the  shock  wave.  The  parameters  of 
the  simulations  are  the  mean  Mach  number  Mj^,  fluctuation  Mach  number  Mt  = 
q/c,  and  turbulent  Reynolds  number  Rca-  described  in  Appendix  resolution 
requirements  limit  the  range  of  the  parameters  in  the  simulation. 

Resolution  of  the  shock  wave  structure  limits  the  range  of  shock  wave  strengths; 
mean  upstream  Mach  number  in  this  work  was  in  the  range  1.05  <  Mf  <  1.20. 

Upstream  of  the  shock  wave,  turbulence  is  weakly  compressible  and  isotropic 
with  Mt  <  0.2,  where  compressibility  effects  are  negligible  [Lee  et  al.  1991b].  The 
range  of  Mt  studied  in  this  work  is  0.057  <  Mt  <  0.110. 

The  resolution  requirement  of  turbulence  length  scales  limits  the  range  of  tur¬ 
bulent  Reynolds  numbers.  Here  we  define  two  turbulent  Reynolds  numbers, 
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where  u'  =  (u  >'/3)>/2  and  the  rate  of  dissipation  of  turbulent  kinetic  energy  per 
unit  volume,  e,  and  the  longitudinal  Taylor  microscale  are  defined  as 
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6.3 


respectively.  In  i!K()inpressi})le  isotropic  turbulence,  the  two  Reynolds  nuinljers 
are  related  as  Rcj  —  The  range  of  the  turbulent  Reynolds  nundjer  in  the 

simulation  was  12  y.  Rc \  25.  corresponding  to  <S0  y-  Rcj  y.  300. 

riius,  we  have  investigated  interactions  between  weak  shock  waves  and  weakly 
compressible  isotropic  turbulence  at  low  Reynolds  numbers.  Table  4.1  lists  the 
simulation  parameters.  (The  values  of  Mf.Rex,  and  Re\  are  taken  at  the  location 
just  l)efore  the  shock.)  In  this  cha])ter,  we  discuss  modification  of  turbulence  by  the 
shock  wave  (Section  4.1 )  and  modificaticjii  of  the  shock  wave  by  turbulence  (Section 
4.2). 


Fable  4.1  Parameters  ior  the  simulations  of  shock-turbulence  interaction 
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4.1  Modification  of  Turbulence 

4.1.1  Preliminary  Considerations 

Figures  4.1(a-c)  show  the  evolutions  of  mean  How  quantities  across  the  shock 
wave.  The  mean  flow  quantities  uiulergo  rapid  jumps  through  the  shock  wave. 
The  average  downstream  values  of  pressure  and  temperature  are  slightly  higher 
tlian  the  values  tor  the  corresponding  laminar  shock  wave,  lliesf  higher  mean 
values  are  cause<l  l)y  the  irreversible  energy  transfer  from  turbulent  kinetic  energy 
to  tlie  internal  energy  lyv  visccuis  di.ssipat ion.  Mean  presstire  and  temperature 


ululergo  slight  overshoots  just  hehiiul  the  shock  wave  tollowecl  l>y  relaxations  due 
to  nonequipartition  of  energy  between  fluctuations  in  pressure  and  velocity,  as  flow 
passes  through  the  shock  wave  [Sarkar  it  al.  1991].  Suppose  the  total  enthalpy  of 
the  flow  «*«*  2  is  constant  behind  the  shock  wave,  then  its  fluctuation 

vanishes,  so  that 

/4'  c]r'  +  -  0. 

The  etteci  of  temperature  overshoot  T'  on  the  velocity  fluctuation  can  be  rep¬ 
resented  as 
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Even  though  the  magnitude  of  the  overshoot  is  small  (about  5%)  compared  to  the 
jump  across  the  shock  wave  with  d/[  -  1.20,  it  contributes  substantially  to  the 

level  of  the  velocity  fluctuations  (about  20%  for  a  flow  with  Mj  —  0.15). 

\elocitv  derivative  skewness  factor  .%»  is  defined  as 
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(4.1) 


which  has  a  value  ol  about  0.4  to  0.0  in  isotropic  turbulence.  This  value  varies 
with  compressiluht y  as  well  as  with  the  turbulent  Reynolds  number  iTavoularis  it 
III.  107-'^.  Erleliaclu’r  it  al.  1990.  Figure  •1.2  shows  tlie  ov(  tlution  of  the  vehx'itv 
derivative  skewnesses  for  ca.se  Turbidence  at  kn-’  \  12  may  be  regarde<l  as 

realistic.  Since  the  mean  position  of  the  shock  wave  is  at  I.S.S,  turbidence 

interacting  with  the  shock  wave  is  considere<l  to  l>e  la-alistic. 

( )ne  dimei\sional  j)ower  s))ect  rum  f.  ^(  i'2- i  )  of  alluctuation  .r  j .  .r.i ,  .  M  from 

the  average  /(j’l  I  is  deflued  by 
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vvluTt'  (■)  denotes  aver<igiii,L:  over  th<'  J";)  dinTtioii  as  well  as  in  time,  /  is  the  P\)urier 
transform  of  /'  in  tliexj  direct  ion,  ami /*  is  its  complex  conjugate.  Figure  4.3  shows 
the  tim'-dimensional  power  spectra  of  velocity  comjionents  and  diuisity  ujjstream 
of  the  shock  wave  for  case  (',  wlu're  k-o/k-c  is  the  largest  {k\-  is  the  cutofT,  or  the 
largest  wave  numher  represented  in  the  simulation).  'I'he  Kolmogorov  wave  number 
k' K  (<  '7'*^’^)'  '  ‘’f  simulation  was  kj^/ko  ~  4.63,  or  kj^/kc  =  0.869.  So,  the 
Kolmogorov  scale  is  actually  capturetl.  Fhe  spectra  decay  at  least  three  orders  of 
magnitude,  which  shows  that  the  How  field  is  well  re.solved.  d’he  spectra  of  £"1  and 
£3  are  in  good  agreement  as  expected  of  isotropic  turbuh'iice.  The  relation  between 
the  .sp('ctra  I''\(k-2)  ami  /''_>(  v-j)  isotropic  turinilence. 


llinze  1!)75  is  also  satisfied. 

Figures  4.4(a-c)  show  the  evolution  of  the  one-dim<’nsional  power  spectra  of 
1/ p  (/!,  and  density,  respi'ct  i  vely,  throughout  the  computational  domain.  Across  the 
shock  wave,  enhancement  of  the  spi'ctra  can  be  noticml  with  more  amplification 
at  large  wave  numbers.  .As  the  How  evolves  furtln'r  downstri'am  of  the  siiock,  the 
spectra  drop  over  tin'  entire  range  of  the  wave  numbers.  More'  amjililication  at 
large  wave  numbers  across  the  shock  wave  leads  to  the  decnuise  of  the  turbulence 
length  scales,  especially  tlu'  4'aylor  microscale. 

I  vvo-point  correlation  and  the  corresponding  integral  scale,  Ayy,^  are 

defined  as 


/'(x)/'(x  t  re, 
/'(x)/'(x) 


‘V/.u  / 

V  0 

respec  I  i  V  el  \ ,  where  the  direction  of  separation  is  indicated  by  o  3,3.  In  order  to 
check  the  ;ide(|uai  V  of  the  computational  liox  size  in  the  transverse  diri-ctions  where 
pf-riodii  lioi.ndaiv  conditions  are  used,  w<'  examine  the  two  point  correlations  of 


three  velocity  components,  density  and  pressure  in  the  X2  direction  upstream  of 
the  shock  for  case  A,  where  ko/kc  is  the  smallest  (Figure  4.5).  The  longitudinzd 
velocity  correlation  C?22,2  decays  monotonically  to  zero  as  expected.  Even  with  some 
problems  in  the  sample  size,  the  lateral  velocity  correlations,  Q\i  2  and  Q33^2i  show 
approximate  isotropy  especiaJly  for  small  separations.  (Correlations  are  calculated 
by  averaging  over  65  saved  fields  which  are  separated  in  time  by  Atu' /Xi  =  0.064.) 
The  correlations  of  pressure  and  density  fluctuations  are  found  to  be  identical.  The 
size  of  the  computationaJ  box  in  the  transverse  directions  appears  to  be  adequate. 

4.1.2  Evolution  of  Turbulent  Kinetic  Energy 

Figure  4.6  shows  the  evolution  of  the  normal  components  of  the  Reynolds  stress 
tensor  72, j  defined  by  [Favre  1965a,  1965b] 


pu'lu'j 


(4.8) 


The  off-diagonal  components  of  R^j  stay  close  to  zero  over  the  entire  flow  field 
due  to  symmetry  (isotropy  upstream  and  axisymmetry  downstream)  in  the  velocity 
fluctuations.  The  statistics  of  the  s*»-eamwise  component  in  the  shock  zone  contains 
the  intermittency  effects  due  to  the  oscillations  of  the  shock.  (The  boundaries  of 
the  shock  oscillations  are  defined  as  the  locations  where  tJi  ]  =  0.)  The  effect  of 
kinematic  oscillation  of  the  plane  shock  wave  on  the  statistics  is  investigated  in 
Appendix  F.  All  the  velocity  fluctuations  are  enhanced  during  the  interaction  as 
predicted  by  the  linear  analyses  (RDT  and  LIA).  Turbulent  velocity  fluctuations  are 
anisotropic  behind  the  shock  wave.  The  return  to  isotropy  is  found  to  be  negligible 
compared  to  the  decay  rate.  The  amplification  in  the  variance  of  transverse  velocity 
fluctuations,  which  is  defined  as  the  ratio  of  the  downstream  maximum  value  to 
the  upstream  minimum  value,  lies  between  the  near-field  and  far-field  predictions 
of  LIA  due  to  the  osci'iatory  movement  of  the  distorted  shock  front:  for  =1-2 
(case  C),  the  simulated  amplification  is  1.19,  while  the  near-field  LIA  prediction 
is  1.45  and  the  far-field  LIA  prediction  is  1.15.  The  streamw’ise  component  /2ii 
undergoes  rapid  increase  behind  the  shock  wave:  the  linear  amdysis  (LIA)  predicts 
monotonic  decay  for  all  the  velocity  fluctuations  (Figures  2. 5-2. 7).  (Discussions 
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on  the  anomalous  behavior  of  the  statistics  iiear  the  outflow  boundary  in  Figures 
4. 6-4. 9,  4.11,  and  4.23  are  given  at  the  end  of  this  section.) 

In  order  to  identify  the  mechanisms  of  amplification  and  rapid  evolution  of  tur¬ 
bulent  kinetic  energy,  the  terms  in  the  transport  equation  of  the  Reynolds  stress 
tensor  R^J  were  computed.  The  transport  equation  for  Rij  [Dussauge  et  al.  1987] 
is 


=  -p  Rik 


du 


dx 


]  n 


(II) 


dx. 


^  dxi 


(III) 


dx , 


^  dxi 
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(4.9) 


The  convection  (I)  of  the  Reynolds  stress  tensor  balances  with  the  production  by 
the  mean  .  ‘rain  field  (II)  and  the  production  by  the  mass  fluctuation  (III),  the 
pressure  work  (IV),  the  turbulent  transport  (V),  and  the  viscous  dissipation  and 
transport  (VT). 

Figure  4.7  shows  the  TKE  budget.  The  statistics  of  the  flow  variables  inside  the 
shock  wave  are  contaminated  by  the  intermittency  effect  caused  by  the  unsteady 
distortion  of  the  shock  wave  (Appendix  F).  Turbulence  amplification  mechanisms 
in  shock-turbulence  interaction  cannot,  therefore,  be  unambiguously  identified  by 
investigating  the  statistics  of  the  numerical  simulations  inside  the  shock  zone.  Out¬ 
side  the  shock  wave,  the  viscous  di.ssipation  is  the  dominant  term  and  the  pressure 
work  term  just  downstream  of  the  shock  is  the  only  other  term  which  has  a  compa¬ 
rable  magnitude.  The  rapid  evolution  of  TKE  just  downstream  of  the  shock  wave 
is  due  to  pressure  work. 

The  pressure  work  term  which  causes  the  rapid  evolution  of  TKE  downstream 
of  the  shock  can  be  decomposed  into  two  terms,  the  pressure  dilatation  p'u" .  and 

the  pressure  transfjort  term  (p'u(*)  ,  as 

7,u;'-=p^.-(p'u»),,.  (4.10) 
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Positive  pressure  dilatation  leads  to  reversible  energy  transfer  from  the  mean  inter- 
naJ  energy  to  the  turbulent  kinetic  energy,  while  the  pressure  transport  is  respon¬ 
sible  for  redistribution  of  TKE  in  the  inhomogeneous  direction.  These  decomposed 
terms  are  shown  in  Figures  4.8(a)  and  (b).  The  decomposition  inside  the  shock 
wave  is  very  similar  to  that  in  the  kinematic  oscillation  of  a  plane  shock  wave  (Ap¬ 
pendix  F),  which  suggests  that  the  behavior  of  the  profiles  inside  the  shock  wave 
are  mainly  due  to  its  unsteady  motion.  The  decomposition  downstream  of  the 
shock  wave,  Figure  4.8(b),  shows  that  the  rapid  evolution  of  TKE  is  caused  by  the 
pressure  transport  term  and  that  the  pressure  dilatation  acts  mainly  to  convert  the 
mean  internal  energy  into  turbulent  kinetic  energy. 

Figures  4.9(a)  and  (b)  show  the  budgets  of  Rn  and  R22  outside  the  shock  wave. 
The  effect  of  the  pressure  work  term  is  quite  pronounced  in  the  equation.  The 
R22  equation  has  no  pressure  transport,  since  the  flow  is  homogeneous  (periodic 
numerically)  in  the  X2  direction.  Therefore,  the  rapid  evolution  of  TKE  is  meiinly 
from  the  streamwise  fluctuations. 

In  numerical  simulations  of  two-dimensional  inviscid  turbulence  interacting  with 
a  shock,  Rotman  [1991]  found  that  turbulence  is  less  amplified  for  the  higher  up¬ 
stream  turbulence  level  for  the  same  mean  shock  strength.  Here,  comparison  is 
made  for  the  amplification  ratios  of  transverse  velocity  fluctuations  of  cases  B  and 
D,  both  with  the  same  shock  strength,  =  1.20.  Turbulent  kinetic  energy  in 
case  B  is  three  times  higher  than  that  in  case  D,  while  turbulent  Reynolds  numbers 
are  comparable.  The  amplification  ratio  is  higher  for  the  weaker  upstream  turbu¬ 
lence  by  about  8%:  the  amplification  ratio  is  1.19  for  case  B  and  1.28  for  case  D. 
Even  though  Rotman’s  simulation  had  a  stronger  shock  wave  with  =  2.07,  the 
reduction  in  amplification  was  about  the  same  for  the  same  change  in  the  upstream 
turbulence  intensity. 

Application  of  Thompson’s  [1987]  boundary  condition  at  the  outflow  generated 
anomalous  statistics  of  the  streamwise  velocity  fluctuations  in  a  region  near  the 
outflow  (Figure  4.6),  pressure  work  (Figures  4. 7-4. 9)  and  dilatation  (Figure  4.11). 
These  anomalies  are  due  to  an  incomplete  suppression  of  the  acoustic  wave  reflec¬ 
tions.  In  order  to  investigate  whether  these  anomalous  behaviors  affect  the  overall 
ev...iution  of  the  flow  downstream  of  the  shock  wave,  the  more  refined  boundary 
conditions  of  Giles  [1990]  was  implemented  in  the  code.  These  numerical  experi¬ 
ments  %^erified  that  these  undesirable  behaviors  were  confined  to  only  a  small  region 
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near  the  boundary.  (The  Giles  boundary  condition  implemented  is  a  first-order  ap¬ 
proximation  for  the  free  propagation  of  a  wave  through  the  boundary,  whereas  the 
Thompson  boundary  condition  is  a  zeroth-order  approximation.)  A  brief  descrip¬ 
tion  of  the  Giles’  boundary  condition  and  a  comparison  of  statistics  are  found  in 
Appendix  G. 

4.1.3  Dilatation  and  Vorticity 

Up  to  the  turbulence  Mach  number  Aft  of  about  0.5,  the  dissipation  rate  of 
turbulent  kinetic  energy  t  can  be  approximated  by  neglecting  the  effect  of  viscosity 
variations,  as  [Lee,  Lele,  and  Moin  1991b] 

=  (^i  +  ,  («i) 

where,  0  =  u^  ^  is  the  dilatation  and  is  the  ratio  of  the  dilatation  and  vorticity 
variances  defined  by 


d'^ 


J.J. 

It 


(4.12) 


Therefore,  the  ratio  is  a  relative  magnitude  of  the  compressible  dissipation  rate 
to  the  incompressible  dissipation  rate. 

The  viscous  term  in  the  TKE  transport  equation,  reduces  to  the  ex¬ 

pression  (4.11)  in  homogeneous  turbulence  if  property  variations  are  neglected. 
This  approximation  is  tested  in  Figure  4.10.  The  contributions  of  turbulence  inho¬ 
mogeneity  and  property  variations  to  the  viscous  dissipation  are  found  negligible 
outside  the  shock  zone. 

As  noted  in  (4.11),  variances  of  dilatation  and  vorticity  contribute  to  the  dissi¬ 
pation  rate  of  turbulent  kinetic  energy.  In  the  following,  we  investigate  the  effect  of 
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a  shock  wave  on  the  variances  of  dilatation  and  vorticity.  The  evolutions  of  dilata¬ 
tion  variance  and  enstrophy  are  shown  in  Figure  4.11.  The  variance  of  fluctuating 
dilatation  is  enhanced  by  three  orders  of  magnitude  across  the  shock  wave,  and 
decays  very  rapidly  behind  the  shock  wave.  Enstrohphy  is  also  amplified  is  passage 
through  the  shock. 

Figure  4.12  shows  the  evolution  of  the  components  of  vorticity  across  the  shock 
wave  for  case  C.  Linear  analyses  predict  that  transverse  components  of  vorticity  are 
amplified  whereas  the  streamwise  component  is  unchanged.  The  simulation  results 
shown  in  Figure  4.12  (and  all  other  simulations)  are  consistent  with  the  linear 
analysis  prediction.  Turbulence  behind  the  shock  wave  becomes  axisymmetric  in 
vorticity  fluctuations  as  well  as  in  velocity  fluctuations. 

Figure  4.13  compares  the  amplifications  of  the  transverse  vorticity  for  different 
shock  wave  strengths.  The  amplification  is  smaller  for  the  weaker  shock.  The  am¬ 
plifications  of  the  transverse  vorticity  across  the  shock  wave  computed  from  the 
simulations  compare  favorably  with  the  predictions  of  linear  analyses:  the  max¬ 
imum  difference  is  5%  for  case  F  with  =  1.05  and  Mt  =  0.10,  where  the 
local  shock  wave  structure  is  significantly  modified  (see  Sec.  4.2.2).  The  effect  of 
turbulent  Reynolds  number  on  the  vorticity  amplification  is  found  to  be  negligible. 

Figures  4.14  (a)  and  (b)  show  the  evolutions  of  vorticity  for  flows  with  different 
turbulent  Reynolds  numbers  interacting  with  a  shock  wave  of  the  same  strength. 
Vorticity  stays  fairly  constant  upstream  of  the  shock,  and  after  the  interaction 
transverse  components  decay  whereas  the  streamw'ise  component  increases  near  the 
shock  wave.  In  order  to  check  if  the  increase  of  the  streamwise  vorticity  component 
behind  the  shock  (Figure  4.14  (a))  is  caused  by  poor  resolution  of  the  simulation, 
a  coarse  grid  simulation  with  97  x  48  x  48  points  was  performed.  Figure  4.14 
(c)  compares  the  statistics  of  vorticity  fluctuations  with  those  from  the  original 
simulation  with  129  x  64  x  64  grid  points.  The  variance  of  the  streamwise  vorticity 
fluctuaticjii  predicted  by  the  coarse  grid  sirnulation  is  slightly  less  than  that  of 
the  original  simulation,  thus,  confirming  that  the  increase  of  the  variance  of  the 
streamwise  vorticity  behind  the  shock  wave  is  not  from  a  numerical  artifact. 

In  order  to  identify  the  dmninant  mechanisms  for  the  vorticity  amplification  and 
the  Reynolds  rmmber  effects  on  the  evolution  of  vorticity.  the  trans])ort  equations 
for  vorticity  variances  were  examined.  The  fluctuating  vorticity  ecjuation,  which  is 
obtained  by  taking  the  curl  of  />  *  ■  (3.11)  is: 
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(4.13) 
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where  s^j  =  2(^1^  +  ^-j^i)  is  the  strain  rate.  Multiplying  (4.13)  by  2u;^  and  taking 
the  average  over  the  homogeneous  transverse  directions  and  over  time  gives: 
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+  '^^a]k^'aP,}P,klp^  -i^a^'a^'k),k  +  (4-14) 
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[Blaisdell  et  al.  1990],  where  the  repeated  Greek  indices  are  not  summed.  Here  4>q 
is  the  viscous  dissipation  and  transport,  given  by 


(4.15) 


The  term  on  the  left  hand  side  of  (4.14)  represents  the  advection  by  the  mean  flow. 
The  first  and  second  terms  on  the  right  hand  side  represent  vortex  stretching  by 
the  mean  and  turbulent  strain  fields,  respectively.  The  next  two  terms  represent 
production  (removal)  by  dilatation.  The  fifth  is  the  vorticity  production  by  the 
baroclinic  torque,  the  sixth  is  the  transport  by  the  turbulent  velocity  field,  and  the 
last  term  is  the  viscous  dissipation  and  transport. 

The  balance  of  the  terms  in  (4.14)  for  the  transverse  vorticity,  u;!,  is  shown 
in  Figure  4.1.5.  The  averages  are  taken  over  the  homogeneous  directions,  X2  and 
X3,  and  over  the  65  saved  spatial  fields  for  case  A  which  are  sepiarated  in  time  by 
Afi//A,  -  (1.064.  All  the  t  erms  in  (4.14)  were  evaluated  independently,  and  the 
maximum  imbalance,  the  difference  between  the  LHS  and  the  RllS,  occurs  outside 
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the  shock  zone  and  is  about  10%  of  the  largest  terms,  the  viscous  dissipation  and 
turbulent  vortex  stretching.  Inside  the  shock  wave,  the  vorticity-dilatation  (IV+V) 
is  the  dominant  source  for  vorticity  amplification.  The  viscous  term  (VIII)  is  bal¬ 
anced  with  the  vortex  stretching  (II+III)  beyond  the  shocked  region.  Baroclinic 
torque  (VI)  is  less  than  1%  of  the  leading  terms  throughout  the  domain  in  all  the 
simulations,  including  case  F  where  the  shock  wave  is  strongly  distorted.  The  effect 
of  turbulent  transport  (VII)  is  also  found  to  be  negligible. 

Figure  4.16  shows  the  vortex  stretching  and  vorticity-compression  by  the  mean 
and  turbulent  velocity  fields,  as  they  appear  in  (4.14).  Even  with  the  overestima¬ 
tion  of  turbulent  strain  rate  in  the  shock  wave  due  to  intermittency,  the  vorticity- 
mean  compression  (IV)  is  much  larger  than  the  other  terms  inside  the  shock.  The 
dominance  of  the  vorticity-mean  compression  explains  the  good  comparison  of  the 
simulation  results  with  the  predictions  of  the  linear  analyses. 

Even  though  amplifications  of  the  transverse  components  of  vorticity  obtained 
in  the  simulations  are  very  close  to  the  predictions  of  linear  analyses,  they  are 
systematically  lower  than  the  linear  predictions.  The  difference  becomes  larger 
for  the  higher  upstream  turbulence  levels,  or  for  the  larger  values  of  the  ratio, 
—  1).  Figures  4.17  (a)  and  (b)  show  the  balance  between  the  nonlinear 

terms  in  the  equation  and  their  net  effect  on  the  amplification.  The  net 
effect  of  the  nonlinear  terms  inside  the  shock  wave  acts  against  the  amplification 
of  vorticity.  However,  this  nonlinear  effects  are  negligible  compared  to  the  linear 
effects  (Figure  4.1.5). 

Figure  4.18  shows  the  balance  of  the  terms  in  (4.14)  for  the  streamwise  vorticity 
■  Inside  the  shock  wave,  effects  of  vortex  stretching  (11  I  III)  and  vorticity- 
dilatation  (IV-fV)  tend  to  cancel  each  other,  resulting  in  no  appreciable  change  in 
the  streamwise  vorticity. 

As  shown  in  F  igures  4.12,  4.14(a)  and  (b),  there  are  differences  in  the  evolution  of 
the  streamwise  vorticity  for  flows  with  different  Reynolds  numbers.  To  identify  the 
effect  of  the  turbulent  Reynolds  number  mi  the  evolution  of  vorticity,  the  vorticity 
budgets  for  flows  with  different  Reynolds  number  are  shown  in  Figures  4.19  (a-d). 
Outside  the  interaction  zone,  the  dominating  terms  in  (4.14)  are  the  viscous  dissi¬ 
pation  (Vni)  and  the  vortex  stretching  (fI  +  111)  maiidy  by  the  turbulent  strain 
rate.  Roth  the  stretching  and  dissipation  increase  by  aliout  the  same  amount  with 


the  increase  of  the  Reynolds  miinher  in  the  budget  of  the  transverse  vorticity  com¬ 
ponent,  resulting  in  no  significant  Reynolds- number  dependence  in  the  evolution 
of  the  transverse  component.  However,  some  Reynolds-number  dependence  was 
found  in  the  budget  of  tlie  -.treamwise  vorticity  evolution:  the  turbulent  stretching 
overtakes  the  viscous  dissipation  behind  the  shock  wave  for  the  higher  Reynolds 
number  flow,  resulting  in  a  higher  increase  in  the  streamwise  vorticity  component 
behind  the  shock  wave  (Figure  4.11(a)). 

I’or  hfjinogeneous  turbulence  with  uniform  density,  the  viscous  dissipation  and 
transport  term  4>a  reduces  to  the  homogeneous  viscous  dissipation 

■S’.','  --  (‘‘•'8) 

I'he  evolutions  of  and  are  nr<'vid<>d  in  Figure  4.20.  The  effects  of  turbulence 
inhomogeneity  and  density  variation  on  tlie  viscous  term  are  found  to  be  negligible 
except  for  4'_)  in  tlu'  sliock  rone. 

4.1.4  Turbulence  Length  Scales 

I'igure  4.21  shows  the  evolution  of  the  integral  scales  defined  in  (4.7)  —  Aj]  2»  ^pp,2i 
and  ■\pp,2-  Across  the  shock  wave,  tl.e  decrease  in  the  transverse  velocity  integral 
scale  A22.2  fpiilt'  pronounced,  d  he  length  scales  of  density  and  pressure  fluctua¬ 
tions  stay  very  close  to  each  other  througliout  the  domain.  They  grow  upstream  of 
the  shcjck  wave,  decrease  .significantly  across  the  shock,  and  recover  to  the  upstream 
levels  in  a  short  distance  downstreani  of  tlie  shock  w'ave. 

fdgure  4.22  shows  the  evfjiution  of  the  turbulence  length  scale  /  defined  as 


.•\s  the  flow  approaches  the  ^hork  wav.*,  this  length  scale  decreases,  probably  because 
the  energy  spec*  ruin  becomes  ftdh  r  at  high  wave  numbers.  (The  inflow  spectrum 
given  in  (3. .42  )  has  virtu, dly  no  energy  b<*yoiid  k / k'o  >  2.)  The  length  scale  decrea.ses 
further  across,  th**  shock  vv.ava  .  an.;  ri  o*;-  re.ther  raj)idly  over  a  short  distance  behind 
the  shock  wace  1  he  lengt  :i  *a,  ai,  -  -t  imat ui  inside  the  shock  wave  is  significantly 
cont  ariuna.t  e,|  b  v  *  ii  ■  over e  -  t  1 1 ; i.a'  n  .11  ■  .f  u'T  j  due  to  t  he  oscillat  ion  of  the  shock  wave. 


Another  length  scale  of  interest  is  the  longitudinal  Taylor  inicroscaJe  defined 
in  (4.2).  Figure  4.23  shows  the  evolution  of  the  Taylor  microscales  throughout  the 
computational  domain.  The  evaluation  of  the  streamwise  microscale  Aj  inside  the 
shock  wave  is  also  significantly  contaminated  due  to  the  shock  wave  oscillation. 
Noticeable  reductions  of  all  the  microscales  are  found  across  the  shock  wave. 

Figure  4.24  compares  the  Taylor  microscale  reductions  in  a  transverse  direction 
from  different  cases  with  predictions  by  LIA  and  RDT.  For  the  range  of  upstream 
Mach  numbers  simulated  =  1.05, 1.IO,  1.20),  the  predictions  by  LIA  and  RDT 
are  virtually  the  same.  Except  for  cases  with  strong  upstream  turbulence  intensity, 
the  simulation  results  compare  favorably  with  the  predictions  by  the  linear  theories. 
As  the  upstream  turbulence  intensity  increases,  the  transverse  Taylor  microsccdes 
decrease  more  across  the  shock  wave  with  the  same  Af^  (compare  cases  A,  B,  and 
C).  The  streamwise  Taylor  microscale  evolves  rapidly  downstream  of  the  shock  and 
it  is  difficult  to  identify  its  “downstream  value.” 

4.1.5  Thermodynamic  Properties 


In  Figure  4.25(a),  we  present  evolutions  of  rms  pressure,  density,  and  tempera¬ 
ture  fluctuations  ipr,Pr,  and  TV).  As  the  flow  passes  through  the  shock  wave,  all 
the  fluctuations  are  amplified,  followed  by  a  rapid  decay.  These  fluctuations  are 
virtually  isentropic.  The  polytropic  exponent  ni(2:i)  is  defined  as 


/  X  Pr(xi  /p(j^l) 

ni(xi)  =  — - — - -, 

pr(xi)/p(xi) 


(4.18) 


which  is  equal  to  7  (1.40  here)  for  isentropic  fluctuations.  Figure  4.25(b)  shows 
evolutions  of  the  polytropic  exponent.  The  polytropic  exponent  stays  close  to  the 
isentropic  value  throughout  the  domain,  varying  between  1.35  and  1.40. 

In  order  to  identify  the  mechanisms  of  amplification  and  decay  of  the  density 
fluctuation  variance  p'  ,  the  budget  for  the  density  fluctuation  variance  [Taulbee  et 
al.  1991]  is  investigated: 
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The  convection  (I)  is  balanced  by  the  productions  by  the  mean  compression  (II) 
and  tlie  mean  density  grailient  (III),  density-dilatation  correlations  (IV-fV),  and 
turbulence  transport  (VI).  Figures  4.26  (a)  and  (b)  show  the  balance  of  the  terms 
in  (4.19)  for  the  density  fluctuation  variance,  fy'^ .  The  averages  are  taken  over  the 
homogeneous  directions,  j-o  and  ^-3,  and  over  the  65  saved  spatial  fields  which  are 
separated  in  time  by  Atu^/A]  —  0.064  for  case  A.  Across  the  shock  wave,  density 
fluctuation.^  are  enhanced  mainly  by  the  production  due  to  the  mean  compression 
(II)  and  the  mean  density  gradient  (III).  Density-fluctuating  dilatation  correlations 
(IV +  \'^)  are  significantly  overestimated  inside  the  shock  wave  due  to  the  shock  front 
oscillation,  but  their  net  effect  is  the  suppression  of  the  density  fluctuation  during 
the  interaction.  Behind  the  shock  wave,  however,  the  evolution  of  is  dominated 
by  the  density-fluctuating  dilatation  correlation  (IV). 

Figure  4.27  shows  the  joint  probability  density  function  of  the  instantaneous 
pressure  versus  instantaueou.s  density  scaled  with  their  local  mean  values,  p{x\) 
and  f){T\  ).  It  is  clear  that  the  isentropic  relations  is  satisfied  for  the  instantaneous 
flow,  even  inside  the  shock  wave.  The  local  polytropic  exponent  n2(x)  is  defined  as 


^2(x) 


p'(x)/p(a-i) 

/’'(x)/p(a'i)‘ 


The  average  polytropic  ex[)onent  rF;>  is  obtained  by  averaging  n2(x)  over  the  flow 
field.  (I  he  states  with  l//(x)/p(  Jj  )i  <  10  ^  were  excluded  in  the  averaging  process 
to  avoid  large  scatter.)  The  average  polytropic  exponent  was  found  to  be  very  close 
to  the  isentropic  value  [n-yh  l.Cl). 

F’.xamining  a  limited  experimental  data  sets,  Morkovin  [1962]  (see  also  Bradshaw 
1977  )  pointed  out  that  in  non-liypersonic  bouiulary  layers  the  acoustic  mode  is 
negligible-  and  the  eutroj'v  mod('  i.s  very  small  for  normal  rates  of  heat  transfer.  lie 
theii  deduced  that 

1.  1 

so  that 
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c)  ii  t 
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(4.20) 


where  the  total  temperature  is  Tt  =  T  +  UiU^/2cp.  This  is  known  as  Morkovin’s 
Strong  Reynolds  Analogy  (the  assumption  of  negligible  total  temperature  fluctu¬ 
ations)  and  is  widely  used  to  correlate  thermodynamic  property  fluctuations  with 
velocity  fluctuations  in  compressible  turbulence  closures.  A  more  general  correla¬ 
tion  through  the  polytropic  coefficient  n  is  suggested  by  Rubesin  [1976]  as 


J  _  rrtit 

P  p  71  I 

p  p  n  —  1  7* 


(4.21) 


Morkovin’s  relation  is  a  special  case  of  (4.21)  with  n  =  0. 

Figure  4.28  shows  the  correlation  between  density  and  temperature  fluctuations, 
Cpj',  defined  by 


CpT  = 


p>T"  p 


which  would  be  n  —  1  if  (4.21)  were  correct.  Since  the  polytropic  coefficient  n 
was  found  to  be  close  to  the  isentropic  exponent  (7  =  1.40),  the  validity  of  the 
Morkovin’s  analogy  appears  to  be  questionable  (with  wrong  signs).  The  same 
conclusion  was  reached  by  Blaisdell  et  al.  [1990]  in  their  numerical  simulation  of 
homogeneous  shear  flow.  Using  Gibbs’  equation  and  the  equation  of  state  of  an 
ideal  gas,  the  pressure  fluctuations  can  be  represented  in  terms  of  the  fluctuations 
in  density  and  entropy  as 


(4.22) 


where  s  is  the  dimensionless  entropy  defined  as  a  =  s* / c^.  When  the  temperature 
(or  entropy)  inhomogeneity  in  the  flow  is  not  significant,  or 
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then  (4.22)  reduces  to  the  polytropic  relations  (4.21)  with  n  =  7.  This  inequality 
can  be  translated  into  terms  which  can  be  estimated  [Thompson  1984,  p.l44]  as 


7P  L 

(4.23) 


The  RHS  of  (4.23)  was  found  to  be  less  than  5%  of  the  dilatation  in  the  present 
simulations,  which  verifies  that  the  relations  between  thermodynamic  property  fluc¬ 
tuations  are  nearly  isentropic. 

Expression  (4.21)  reduces  to  the  Strong  Reynolds  Analogy  (4.20)  if  there  exists 
an  appreciable  mean  temperature  gradient  in  the  flow  and  the  pressure  fluctuations 
are  negligible.  The  Strong  Reynolds  Analogy  is,  therefore,  a  good  approximation 
in  turbulent  boundary  layers  where  mean  gradients  of  temperature  and  density 
normal  to  the  wall  are  large. 

4.1.6  Modeling  Issues 

In  the  k-e  formulation  for  compressible  turbulence,  the  transport  equation  for 
turbulent  kinetic  energy  (equation  (4.9)  with  i  =  j)  has  more  rnodeled  terms  than 
its  incompressible  counterpart.  The  additional  terms  are  the  compressible  dissipa¬ 
tion  6^  (included  in  the  pressure-dilatation  correlation  p'u'f-  (included  in 

-u'Ip'^)  and  the  average  turbulent  mass  flux  u" .  Zeman  [1990]  and  Sarkar  et  al. 
[1991]  proposed  models  which  parametrize  compressible  dissipation  in  terms  of  in¬ 
compressible  dissipation  e^,  and  the  fluctuation  Mach  number  Aft  as  f{Mt)- 

Both  models  were  successful  in  predicting  the  suppression  of  the  spreading  rate  of 
compressible  mixing  layers  at  high  convective  Mach  numbers.  Coleman  and  Man- 
sour  [1991]  proposed  a  modified  turbulence  model  for  the  incompressible  dissipation 

when  turbulence  is  subject  to  mean  compression. 

Zeman  [1991a]  has  identified  the  pressure-dilatation  correlation  with  the  rate 
of  change  of  compressible  potential  energy,  represented  by  the  pressure  fluctuation 
variance  p'^ ,  as 


p'Kt 


1  uj  dp'^ 


(4.24a) 
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Zeman  assumed  that  the  rate  of  change  of  is  governed  by  a  nonlinear  relaxation 

A 

mechanism  which  drives  to  an  equipartition  vedue  [Saxkar  et  al.  1991]: 


-  .  -  pI 

'  9x,  T„ 

where  the  acoustic  time  scale  is 


MlT 

+  M2/3) 


(4.246) 


(4.25a) 


(t  =  pfp"  denotes  the  turbulence  time  scale),  and  the  equilibrium  value  Pg  is 


A  _  oM}  +  0Mt  ,,,, 

^  l+aU} +3M>'  '  ' 

a  =  1  and  /?  =  2  were  chosen  to  best  match  the  DNS  results  for  highly  com¬ 
pressible  turbulence  [Plaisdell  et  al.  1990).  The  expression  (4.256)  combines  two 
assumptions:  (1)  the  equilibrium  ratio  of  compressible  to  solenoidal  turbulent  ki¬ 
netic  energy  is 


ql/qj  =  aA/2  + 

and  (2)  in  equilibrium,  the  compressible  potential  and  kinetic  energies  are  in 
equipartition, 


Pc 


-2  2-2 
P  ■ 


Zeman  [1991b]  and  Durbin  and  Zeman  [1991]  proposed  modifications  for  (4.24a) 
and  (4.246)  to  account  for  the  effect  of  the  mean  strain  rate  as 


1  uj  dp'^  p'"  flu, 

2  p  0x|  p  (9x, 


(4.26n) 
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and 
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-  2(7  +  cj)  p‘ 
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dxi' 


(4.266) 


where  =  (5  -  37)/12. 

As  shown  in  Figures  4.29(a)  and  (b),  the  model  in  (4.26a)  appears  to  be  very 
accurate. 

Figure  4.30  shows  an  evaluation  of  (4.266)  using  DNS  data.  The  closure  equa¬ 
tion  for  p'^  was  found  to  be  inadequate  throughout  the  domain.  The  equilibrium 
pressure  variance  Pg  which  best  matches  the  upstream  and  downstream  evolution 
of  was  found  to  be  0.3  times  the  expression  given  in  (4.246),  and  the  result 
is  also  shown  in  the  figure.  Zeman  assumed  that  the  pressure  variance  p'^  is  the 
compressible  potential  energy^  and  relaxes  to  the  equilibrium  compressible  pres¬ 
sure  variance  pf.  Investigation  of  the  simulation  database  of  decaying  compressible 
turbulence  [Lee  et  al.  1991b]  showed  that  the  contribution  of  the  incompressible 
pressure  is  appreciable  even  in  fairly  compressible  turbulence  with  Mt  ~  0.5.  For  a 
flow  with  lower  Mt,  the  incompressible  pressure  comprises  the  major  portion  of  the 
pressure  fluctuations.  Therefore,  the  assumption  behind  the  expression  for  Pg  given 
in  (4.256)  is  invalid  except  for  highly  compressible  turbulence.  For  the  range  of  pa- 
rcuneters  studied  in  this  work,  compressibility  effects  are  not  significant  anywhere 
except  in  the  shock  zone. 

Sarkar  [1991]  has  Jilso  proposed  a  model  for  the  pressure-dilatation  correlation 
using  the  DNS  database  of  isotropic  turbulence  and  sheared  homogeneous  turbu¬ 
lence.  He  developed  his  model  using  the  statistics  in  the  fully  developed  stage 
where  turbulent  kinetic  energy  decays  (isotropic  turbulence)  or  grows  (sheared  tur¬ 
bulence)  in  time.  Taulbee  and  Van  Osdol  [1991]  proposed  a  single  expression  for  the 
sum  of  the  pressure-dilatation  correlation  and  compressible  dissipation,  p'u'f  ■  -I-  c^, 

•1* 

and  tuned  model  parameters  to  match  the  experimental  data  in  a  boundary  layer 
and  a  mixing  layer.  Both  models  were  unable  to  predict  the  present  observation  of 
reversible  energy  transfer  from  mean  thermal  energy  into  turbulent  kinetic  energy 
in  the  absence  of  mean  strain  rate  behind  a  shock. 
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The  pressure  transport  term,  — was  identified  to  be  the  driving  mecha¬ 
nism  for  the  rapid  evolution  of  the  velocity  fluctuations  downstream  of  the  shock 
wave.  Assuming  the  isentropic  relation  between  thermodynamic  fluctuations,  the 
pressure  velocity  correlation  can  be  expressed  as 


—p'u'l  =  —  =  —c^p'u'l. 

1  P  I  I 


(4.27a) 


The  accuracy  of  this  assumption  is  checked  in  Figure  4.31(a),  which  shows  that  he 
use  of  isentropic  thermodynamic  relations  in  (4.27a)  is  satisfactory  throughout  the 
domain.  Therefore,  modeling  of  the  pressure  transport  can  be  reduced  to  that  of 
the  turbulent  mass  flux,  p'u'^.  Taulbee  et  al.  [1991]  developed  a  model  transport 
equation  for  p'u'^  using  Morkovin’s  hypothesis.  However,  as  was  shown  in  section 
4.1.4  applicability  of  the  Morkovin’s  hypothesis  in  the  absence  of  mean  temperature 
gradient  is  questionable.  Zeman  [1991b]  proposed  a  model  rate  equation  for  p'u" 
as 
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(4.276) 


which  drives  the  mass  flux  fluctuation  to  zero  except  near  the  shock,  on  the  fast 
acoustic  time  scale  Ta-  The  accuracy  of  the  model  in  (4.276)  is  checked  in  Figure 
4.31(b)  which  shows  that  it  qualitatively  represents  the  behavior  downstream  of  the 
shock  wave.  However,  there  are  differences  in  the  peak  position  and  its  magnitude. 

To  gain  further  insight  into  the  physics  of  the  turbulent  mass  flux  p'u'^  down¬ 
stream  of  the  shock,  several  different  scalings  were  examined.  The  best  scaling 
is  based  on  the  assumption  that  the  correlation  between  the  density  and  velocity 
fluctuations  is  composed  mciinly  from  acoustic  waves.  The  streamwise  distance  and 
the  turbulent  mciss  flux  are  scaled  as 


ij/ui  p'u'{ 

— — —  and  - , 

A2/C  p  cM^ 


(4.28) 


respectively.  Figure  4.32(a)  and  (b)  show  the  unsealed  and  the  scaled  turbulent 
mass  flux  terms  for  different  upstream  conditions  and  shock  strengths  (M^  — 
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1.10, 1.20),  respectively.  The  match  in  amplitudes  and  length  scales  of  the  turbulent 
mass  fluxes  from  different  simulations  is  significantly  improved  by  the  seeding  (4.28), 
which  suggests  that  the  rapid  evolution  of  TKE  may  have  been  caused  by  the 
propagation  of  acoustic  waves  which  are  generated  during  the  interaction. 


4.2  Modification  of  a  Shock  Wave 
4.2.1  Statistics  of  a  Shock  Wave 

The  characteristics  of  the  shock  wave —  shock  wave  thickness  and  shock  front 
distortion —  vary  with  time  due  to  the  effect  of  upstream  turbulence.  The  peak 
compression  ^min(^2)^3)0  inside  the  shock  wave  along  the  mean  streamlines  is 
used  as  a  measure  of  the  shock  wave  strength.  (The  ^lin  i®  taken  along  the  ij- 
direction  for  each  X2  and  13.  d^iin  is  necessarily  confined  to  a  plane,  however 
for  clarity  its  “stretched”  version  on  a  plane  will  be  presented.)  Figure  4.33  shows  a 
typical  contour  plot  of  ^min(®2) ®3> 0-  (Dashed  contours  denote  values  with  larger 
than  the  corresponding  laminar  value  at  the  same  upstream  Mach  number.)  The 
average  peak  compression  is  found  to  decrease  from  the  peak  compression  of  the 
laminar  shock  wave  by  about  10%.  The  peak  compression  varies  widely  across  the 
transverse  plane,  which  is  reflected  in  the  large  value  of  the  ratio  between  the  rms 
and  the  mean  values  of  the  peak  compression,  (^inin)rms/|^minl  0.42. 

Figures  4.34(a)  and  (b)  show  the  probability  density  function  (PDF)  of  the  peak 
compression  inside  the  shock  wave.  The  flow  tends  to  have  frequent  events  of  large 
compression.  This  trend  is  clearly  shown  in  Fig.  4.34(b),  where  the  probability  of 
large  compression  zones  is  higher  than  the  Gaussian  distribution  by  several  orders 
of  magnitudes.  This  is  confirmed  by  the  skewness  and  flatness  values  of  the  PDF, 
—  0.81  and  11.0,  respectively. 

The  statistics  of  the  shock  front  distortion  were  estimated  by  LIA  in  Section  2.2. 
The  shock  wave  front  in  the  simulated  field,  however,  is  not  clearly  defined,  because 
in  the  numerical  simulations  the  shock  wave  spans  over  several  grid  points.  The 
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pressure  half-rise  point  was  chosen  to  be  the  shock  front  position,  ^{x2,x^,i).  In 
other  words, 


^(i2,i3,0 


Pii,X2,X2,t) 


(4.29) 


where  pjr  denotes  the  pressure  in  the  laminar  shock  wave  which  was  shown  in  Figure 
4.1  to  be  near  the  mean  turbulent  Vcdue.  This  designator  is  very  well-defined  and 
remains  relatively  noise-free  for  the  cases  with  weak  upstream  turbulence.  In  order 
to  check  the  sensitivity  of  the  shock  front  statistics  to  the  special  choice  of  the 
designator,  the  contour  plot  of  the  shock  wave  position  based  on  the  pressure  half¬ 
rise  point  was  compared  with  that  based  on  the  density  half-rise  point  in  Figures 
4.35(a)  and  (b).  The  difference  between  the  rms  shock  distortions  obtained  by  the 
two  methods  is  always  less  than  1%  of  the  predicted  rms  values. 

Figure  4.36  shows  the  seeded  rms  displacement  of  the  shock  front  from  different 
simulations.  The  scaling  suggested  by  LIA  was  found  to  collapse  the  simulation  data 
reasonably  well.  The  LIA  prediction  of  the  shock  front  displacement  is  dependent 
on  the  shape  of  the  spectrum,  especially  on  the  low  wave  number  part.  The  rms 
displacements  of  the  shock  front  from  the  simulation  are  systematically  lower  than 
the  LIA  prediction  based  on  the  energy  spectrum  (2.1)  and  higher  than  that  based 
on  the  von  Karman  energy  spectrum  (2.2).  This  can  be  attributed  to  the  shape  of 
the  upstream  spectrum  in  the  simulation:  As  the  spectrum  develops  from  the  inflow 
spectra  (3.32)  or  (2.1),  it  loses  energy  at  small  wave  numbers  and  gains  energy  at 
large  wave  numbers.  The  scaled  shock  front  displacement  is  smaller  for  the  higher 
upstream  turbulence  level. 

Figure  4.37  shows  the  scaled  rms  shock  front  inclination  angle  in  the  Z2‘<iirection, 
<T2  =  d^/dx2-  LIA  predicts  the  rms  shock  front  inclination  angle  to  be  independent 
of  the  upstream  energy  spectrum  shape.  The  statistics  from  the  simulation  are  in 
fair  agreement  with  the  LIA  predictions.  As  the  fluctuation  Mach  number  of  the 
simulation  increases,  the  simulation  result  deviates  further  from  the  linear  predic¬ 
tion  (Note  that  an  infinitesimal  fluctuation  Mach  number  is  one  of  the  assumptions 
for  valid  application  of  the  linear  anedysis). 

Shock  front  curvature  profiles  obtained  from  the  simulation  was  too  noisy  to 
permit  any  decisive  conclusion  except  for  case  C  where  the  shock  wave  is  resolved 
best.  The  rms  values  of  the  scaled  shock  front  curvature 
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range  from  1.86  to  2.03,  while  the  LI  A  prediction  with  spectrum  (2.1)  is  1.41  for 

Mf  =  1.20. 

4.2.2  Instantaneous  Fields 

Instantaneous  density  fields  at  a  typical  xj  —  X2  plane  from  the  cases  D  and 
F,  are  given  in  Figures  4.38(a)  and  (b):  Figure  4.38(a)  is  for  a  weak  upstream 
turbulence  {M^  =  1.20,  =  0.057)  and  Figure  4.38(b)  is  for  a  relatively  intense 

turbulence  {M^  =  1.05,  =  0.10).  The  overlaid  contour  lines  near  the  center 

of  the  figures  are  iso-compression  (same  Uj  j)  lines.  Figure  4.38(a)  shows  a  clear 
shock  front  across  which  significant  increase  in  density  is  noticed.  For  the  more 
intense  upstream  turbulence  (Figure  4.38(b)),  the  shape  of  the  shock  front  is  more 
distorted.  The  variation  of  peak  compressions  inside  the  shock  wave  in  transverse 
directions  becomes  stronger  for  the  more  intense  upstream  turbulence.  The  varia¬ 
tion  in  the  visual  thickness  of  the  shock  wave  is  also  larger  for  the  stronger  upstream 
turbulence.  Low  density  regions  are  often  found  behind  the  mean  shock  position 
for  flows  with  Mi  >  —  1. 

By  comparing  Figures  4.33  and  4.35(a),  we  found  a  moderate  correlation  between 
the  peak  compression  and  the  shock  front  distortion:  the  peak  compression  is  larger 
for  the  shock  wave  pushed  downstream  relative  to  the  mean  position  by  upstream 
turbulence  and  weaker  for  the  one  pulled  upstream.  Figures  4.39(a)  and  (b)  show 
the  profiles  of  dilatation  along  the  xj-direction  for  the  cases  D  and  F,  respectively. 
The  strength  of  the  shock  wave  (or  the  peak  compression  in  the  shock)  varies  widely 
from  one  streamline  to  another.  For  the  case  of  strong  upstream  turbulence  (Figure 
4.39(b)),  the  structure  of  the  shock  wave  is  significantly  modified:  multiple  peaks 
in  compression  along  streamlines  (e.g.,  o)  are  noticeable.  Each  compression  peak 
has  a  strength  comparable  to  that  of  the  laminar  shock  wave.  Sometimes,  a  shock 
wave  is  replaced  by  a  series  of  compression  waves  (e.g.,  x). 

Figures  4.40(a)  and  (b)  show  the  pressure  profiles  along  the  xj-direction.  For  rel¬ 
atively  weak  upstream  turbulence,  pressure  rises  monotonically  from  the  upstream 
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to  the  downstream  value.  On  the  other  hand,  for  strong  upstream  turbulence,  pat¬ 
terns  of  the  pressure  rise  are  varied:  rapid  monotonic  rise,  slow  monotonic  rise  and 
multiple-staged  rise. 

Figure  4.41  shows  the  streamwise  Mach  number  in  a  near  the  shock 

wave  for  case  A.  The  streamwise  Mach  number  is  a  good  representative  of  the  up¬ 
stream  shock  normal  Mach  number  in  cases  where  <72  <  10°.  The  drops  in  the 
streamwise  Mach  number  across  the  shock  wave  along  the  xj -direction  are  more  or 
less  uniform:  the  higher  upstream  Mach  numbers  correspond  to  the  higher  down¬ 
stream  Mach  numbers,  and  vice  versa.  For  a  shock  wave  fixed  in  space,  however, 
higher  upstream  Mach  number  would  correspond  to  lower  downstream  Mach 
number  M^,  as 


(Mff  = 


1  +  ^(Mf  )2 


(4.30) 


The  uniform  Mach  number  drop  for  different  upstream  Mach  numbers  implies  that 
the  local  fluctuating  shock  wave  speed  tends  to  be  in  phase  with  the  upstream  Mach 
number —  positive  for  higher  Mach  numbers  and  negative  for  lower  Mach  numbers, 
so  that  the  effective  shock-normal  Mach  numbers  are  more  or  less  uniform  across 
the  transverse  plane.  The  nonuniformity  in  the  upstream  Mach  number  is  smoothed 
out  by  the  fact  that  the  local  shock  wave  speed  tends  to  be  in  phase  with  that 
of  the  fluctuating  velocity  (Section  2.2). 
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Figure  4.4.  Evolution  of  the  energy  spectra  of  (a)  streamwise  velocity,  (6)  trans¬ 
verse  velocity,  and  (c)  density  fluctuation  for  case  C: - upstream(/:oX|  =  15.8), 

- downstream(fco3;i  =  21.9), - at  exit  boundary(fcoa^l  =  37.7). 


Figure  4. 


4.  continued. 
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for  case 


Figure  4.7.  Budget  of  terms  (scaled  with  poCo^"o^o)  in  the  turbulent  kinetic 
energy  transport  equation  for  case  C.  (a)  near  the  shock  wave  and  {b)  in  the  entire 
domain:  -a - convection(I), - production  by  the  mean  strain(II), - pro¬ 
duction  by  mass  flux  fluctuation(III),  .  pressure  work(IV),  -  turbulent 

transport(V), - viscous  dissipation  and  transport(VI).  (Vertical  dashed  lines 

denote  the  boundaries  of  shock  intermit tency.) 

94 


ty 

Figure  4.9.  Budget  of  terms  (scaJed  with  po<^o'^o^o)  in  the  transport  equation  of 
Reynolds  stresses,  (a)  Rjj  and  (6)  R221  - convection(I), - pro¬ 
duction  by  the  mean  strain(II), - production  by  ma.ss  flux  fluctuation(III), 

pressure  work(IV), -  turbulent  transport(V), - viscous  dissipation 

and  transport(VI). 
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Figure  4.10.  Approximation  of  the  viscous  term  (scaled  with  po^o'^o^o)  for  case 
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Figure  4.12.  Evolution  of  components  of  vorticy  for  case  C  (Rej^  =  84.8):  - wj, - u;2, 
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Figure  4.15.  Budget  of  terms  (scaled  with  Cou^/fco)  in  the  transport  equation  of  for  case  A:  -e -  con- 

vection(I), - stretching(II+III), - vorticity-dilatation(IV+V),  .  baroclinic  torque(VI), - turbulent 

transport(  VII), - viscous  dissipation  and  transport(VIII). 
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Figure  4.16.  Production  of  vorticity  by  the  mean 
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Figure  4.17.  IndividuaJ  (a)  and  tottd  effects  (b)  of  nonlinear  terms  (scaled 
2  . 2 

with  Couffho)  in  the  transport  equation  of  for  case  A.  In  Figure  4.1 7(o): 

-  stretching(III), - vorticity-compression(V), - baroclinic  torque(V'I), 

transport( VII), - viscous  dissipation  and  transport( VIII). 
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■e - convection(lj, - stretcmng(ll+lllj, - vorticity-compression(lV  +  Vj, .  barochnic  torquei 

turbulent  transport(VII), - viscous  dissipation  and  transport(VIII). 
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Figure  4.19.  Budget  of  terms  (scaled  with  cpu^ko)  in  the  transport  equation 

of  (a)  {Rej'  =  240),  (6)  a;^2  ~  240),  (c)  ~  84.8),  and  (d) 

(jj^2  ~  -  convection(I), - stretching(II  +  III), - vorticity- 

compression(IV+V),  baroclinic  torque(VI), - turbulent  transport(VII), 

- viscous  dissipation  and  transport(VIII). 
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Figure  4.22.  Evolution  of  the  turbulence  length  scale  defined  in  (4.17)  for  case  C.  (Vertical 
boundaries  of  shock  intermittency.) 
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Figure  4.25.  Evolutions  of  (a)  rms  values  of  therodynaniic  properties  and 
(6)  the  polytropic  exponent  rq  for  case  A.  For  (a):  -  P^/Pin’ - P^/Pin’ 
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Figure  4.26.  Budget  of  terms  (scaled  with  {po^ft)‘^^oko)  in  the  transport  equa- 
— 

tion  of  for  case  A:  (a)  in  the  shock  zone,  (6)  in  the  entire  domain,  -b - convec- 

tion(I), - production  by  the  mean  dilatation( II), - production  by  the  mean 

density  gradient(III),  turbulent  density-dilatation  corrclation(IV), - tur¬ 
bulent  density  squared-dilatation  corrclation{ V), - turbulent  transport( VI). 
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Figure  4.28.  Evolution  of  the  correlation  coefficient  c„t  between  density  and  temperature  fluctuati 


koXl 

Figure  4.29.  Evaluation  of  the  Durbin  and  Zeman’s  [1991]  modeling  idea  for 
p^u'f  (a)  near  the  shock  and  (6)  in  the  entire  domain  for  case  C;  -  LHS, 
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- RHS. 
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boundaries  of  shock  intermittency.) 
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Figure  4.31.  Evaluation  of  the  models  (4.27)  for  \  for  case  C:  (a) 

approxmation  in  terms  of  p'u'^  ( -  LHS  of  (4.27a), - RHS  of  (4.27a));  (b) 

transport  of  p'u'^  ( - LHS  of  (4.276), - RHS  of  (4.276)). 
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Figure  4.32.  Unsealed  (a)  and  scaded  (6)  turbulent  mass  fluxes  is  the  position 

for  peak  negative  turbulent  mass  flux):  - case  A, - case  B, - case  C, 

.  case  D, - case  E. 
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laminar  shock  wave,  6  =  —8.1.  Increment  between  lines, 
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Figure  4.34.  Probability  density  function  for  the  peak  negative  dilatation  in 

(a)linear-linear  and  (6)linear-logarithmic  coordinates  for  case  C:  -  PDF  of 

^min’ - PDF  of  a  Gaussian  distribution. 
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upstream. 
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o  case  D,  V  case 
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lines  near  the  center  of  the  figure  are  iso-compression  lines. 
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Figure  4.41.  Contour  plot  for  the  streamwise  Mach  number  for  case  A.  Solid  lines  denote  supersonic  condition, 
and  dashed  lines  denote  subsonic  condition.  Increment  between  lines  AA/i  =  0.05. 
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CHAPTER  5 


CONCLUSIONS  AND  RECOMMENDATIONS 

This  work  has  been  concerned  with  a  numerical  study  of  the  interaction  of  tur¬ 
bulence  with  a  shock  wave.  The  methods  used  were  linear  analyses  (RDT  and  LIA) 
and  direct  numerical  simulation  of  the  compressible  Navi**r  Stokes  equations.  Lin¬ 
ear  analyses  were  used  to  predict  modifications  in  turbulence  statistics  in  passage 
through  the  shock  and  shock  front  statistics.  In  direct  numerical  simulation,  the 
full  equations  were  solved  by  an  explicit  code,  with  the  spatial  derivative  evaluated 
by  a  modified  Fade  scheme.  A  summary  of  conclusions  is  provided  below,  followed 
by  recommendations  for  future  work.  Conclusions  are  divided  into  three  main  ar¬ 
eas:  spatially  decaying  turbulence,  turbulence  modification  by  shock  turbulence 
interaction,  and  shock  wave  modification  by  the  interaction. 


Spatially  Decaying  Turbulence 

A  method  of  generating  stochastic  inflow  boundary  conditions  with  prescribed 
spectrum  was  developed.  Turbulence  intensity,  rms  vorticity,  and  velocity  deriva¬ 
tive  skewness  factor  compare  favorably  with  those  from  the  temporal  simulation. 
However,  the  statistics  of  dilatation  show  significant  departure  from  those  obtained 
in  the  temporal  simulation.  Because  of  this  difference,  caution  must  be  exercised 
in  using  periodic  (or  temporal)  simulation  databases  to  examine  compressibihty- 
driven  quantities  such  as  dilatation  dissipation  and  pressure  dilatation  correlation. 
Turbulence  statistics  from  spatially  evolving  simulations  with  low  cornpressibifity 
effects  are  also  found  to  be  in  agreement  with  the  experimental  data. 


Turbulence  Modification  by  Shock  Turbulence  Interaction 

The  simulations  and  linear  analyses  of  shock  turbulence  interaction  show  that  the 
normal  components  of  the  Reynolds  stresses,  Raa-,  fi^re  enhanced  across  the  shock 
wave.  The  stress  amplifications  are  larger  for  the  stronger  shock  w'aves  within  the 
range  of  Mach  numbers  in  the  simulation,  1.05  <  A/[  <  1.20.  LIA,  however, 


predicts  less  amplification  of  the  streamwise  turbulence  intensity  for  the  stronger 
shock  waves  when  >  2.0,  a  range  beyond  that  of  the  present  simulations. 

The  simulations  show  a  rapid  evolution  of  Raa  immediately  behind  the  shock 
wave,  including  an  increase  in  value  just  downstream.  LIA,  on  the  other  hand, 
predicts  that  Raa  monotonically  decays  from  its  post-shock  value  to  the  far-field 
value.  The  budget  of  the  TKE  transport  equation  revealed  that  the  pressure  work 
term  —u'!p'^  is  responsible  for  this  rapid  evolution.  By  decomposing  the  pressure 

work  term  into  the  pressure-dilatation  correlation  p'u'^ ^  and  the  pressure  transport 

—  {p'u'-)  {  term,  it  was  found  that  the  pressure  transport  is  the  main  contributor  to 
the  pressure  work  term.  Therefore,  rapid  evolution  of  TKE  is  caused  mainly  by  the 
redistribution  of  turbulent  kinetic  energy  in  the  streamv.'ise  direction  at  the  early 
stages  of  relaxation  from  compression. 

The  simulations  and  linear  analyses  show  that  the  rms  of  the  transverse  vortic- 
ity  components  are  amplified  and  that  the  streamwise  vorticity  is  not  influenced 
by  the  interaction.  The  amplifications  are  larger  for  stronger  shock  waves.  The 
amplification  of  the  transverse  vorticity  components  predicted  by  the  simulations 
are  in  excellent  agreement  with  those  of  linear  analyses  for  Mt  <  —  1.  For 

Mt  >  —  1,  the  vorticity  amplification  obtained  from  the  sii  lulation  is  signif¬ 

icantly  less  than  the  estimates  of  the  linear  analyses.  Examination  of  the  budget 
of  u;'^  revealed  that  the  vorticity-compression,  is  the  main  contributor  to 

the  transverse  vorticity  amplification  during  the  interaction.  In  most  cases,  the 
baroclinic  torque  is  at  least  two  orders  of  magnitude  smaller  than  the  vorticity- 
compression  term  in  the  interaction  zone. 

Turbulent  length  scales  decrease  through  the  interaction.  Linear  analyses  predict 
the  amplification  of  one-dimensional  energy  spectra  during  the  interaction,  both  in 
frequency  and  wave  number  space.  In  the  freqiiency  spectrum,  amplification  is 
larger  for  the  smaller  frequencies.  In  the  wave  number  spectrum,  however,  amplifi¬ 
cation  is  larger  at  the  larger  wave  numbers.  More  enhancement  at  the  smaller  scales 
leads  to  a  decrease  in  turbulence  length  scales.  Numerical  simulations  confirmed 
that  the  energy  spectrum  is  amplified  more  at  the  larger  wave  numbers,  whereas 
turbulence  length  scales  are  decreased  through  the  interact if)n. 


We  investigated  the  relations  between  the  fluctuations  in  pressure,  density,  and 
temperature  by  the  use  of  the  simulation  results.  Morkovin’s  hypothesis  was  tested 
and  found  inaccurate.  Instead,  the  relations  between  properly  scaled  rms  property 
fluctuations  are  very  close  to  isentropic  at  least  for  turbulence  passing  through 
a  weak  shock  wave,  <  1.20.  Isentropic  relations  are  satisfied  even  for  the 

instantaneous  fluctuations  throughout  the  flow  field  including  the  shock  wave. 


Shock  Wave  Modification  by  Shock  Turbulence  Interaction 

Due  to  the  nonuniformity  of  upstream  turbulence,  the  shock  wave  has  a  time- 
dependent  distorted  front  and  nonuniform  thickness  in  the  transverse  directions. 
For  the  simulations  with  Mt  <  —  1,  shock  waves  have  well-defined  fronts  with 

single  compression  peaks  in  the  streamwise  direction. 

Through  LI  A,  it  was  found  that  the  shock  front  distortion  scales  with  the  up¬ 
stream  integral  length  scale  and  turbulence  intensity.  We  also  found  that  the  local 
shock  front  inclination  angle  scales  only  with  upstream  turbulence  intensity,  and 
the  shock  front  curvature  is  scaled  with  turbulence  intensity  and  the  inverse  of  the 
Taylor  microscale.  All  scaling  factors  depend  only  on  the  shock  strength.  The 
statistics  of  the  shock  front  obtained  by  the  simulations  are  in  fair  agreement  with 
the  LIA  predictions. 

Through  LIA,  local  shock  wave  speed  was  found  to  scale  with  the  upstream 
fluctuating  velocity.  Instantaneous  flow  fields  from  the  simulations  showed  that  the 
velocity  jump  across  the  shock  wave,  or  the  shock  wave  strength,  is  more  or  less 
uniform  in  the  transverse  directions:  The  shock  wave  moves  upstream  for  lower 
velocity,  and  moves  downstream  for  h'gher  velocity  of  the  approaching  flow,  thus 
reducing  the  variation  in  the  effective  shock-normal  Mach  number. 

For  a  simulation  with  Mt  >  —  1,  shock  waves  no  longer  have  well-defined 

fronts  in  the  transverse  directions.  Low  density  regions  are  often  found  behind  the 
mean  shock  position,  and  shock  wave  thickness  varies  quite  widely  in  transverse 
directions.  Along  the  streamwise  direction,  multiple  peaks  in  dilatation  and  pres¬ 
sure  are  found,  which  are  similar  to  those  observed  in  the  experiments  of  a  very 


weak  shock  wave  interacting  with  a  highly  inhomogeneous  medium  [Hesselink  et  al. 
1988]. 


A  comprehensive  quantitative  data  base  has  been  generated  which  may  be  used 
to  develop  and  test  turbulence  models  and  to  further  study  the  physics  of  shock/tur¬ 
bulence  interaction. 


Recommendations  for  Future  Research 

Direct  numericcil  simulation  of  shock/turbulence  interaction  is  a  new  area  of  re¬ 
search,  and  many  questions  remain  to  be  answered.  From  our  work,  we  recommend 
the  following  directions  for  future  research. 

Interactions  with  Significant  Acoustic  or  Entropy  Fluctuations 

In  this  work,  upstream  turbulence  conditions  in  the  simulations  were  restricted 
to  quasi-incompressible  states.  It  is  also  important  to  study  the  interaction  of  a 
shock  wave  both  with  highly  compressible  turbulence  and  with  flows  of  large  density 
variations. 

Towards  Turbulent  Boundary  Layer  Interaction  with  a  Shock  Wave 

The  logical  next  step  towards  understanding  of  shock/turbulent  boundary  layer 
interaction  is  to  study  interaction  of  homogeneous  turbulence  under  a  mean  shear 
with  a  shock  wave.  Through  this  study,  one  can  investigate  the  effect  of  a  shock 
wave  on  turbulence  anisotropy  and  Reynolds  shear  stress. 

Large-Eddy  Simulation  of  Shock  Turbulence  Interaction 

Direct  numerical  simulation  of  shock  turbulence  interaction  is  restricted  to  low 
Reynolds  numbers  and  weak  shock  waves  by  the  resolution  requirements  of  the 
shock  wave  structure.  Subgrid-scale  models  have  been  actively  under  development 
[Erlebacher  et  al.  1990,  Moin  et  al.  1991]  for  the  large-eddy  simulation  of  com¬ 
pressible  turbulence.  For  shock  turbulence  interaction,  the  ratio  of  the  turbulence 
length  scale  I  —  pq^ / e  to  the  shock  wave  thickness  bg  is 


i:i6 


I  Rc'f 

typically  0(10^)-C1(10^)  or  larger  in  practical  flows  where  Re^p  >  1000.  Since  the 
cost  of  resolving  the  shock  wave  structure  in  large-eddy  simulation  is  prohibitive, 
one  needs  a  subgrid-scale  model  for  the  shock  wave  effect  on  turbulence  [Zeman 
1991b],  as  well  as  for  the  effect  of  small  sccde  turbulence  on  the  large  scales.  Another 
possible  option  is  to  use  a  shock  capturing  technique  [Yee  1987,  Harten  al.  1987] 
without  any  subgrid-sccde  model  of  the  shock  wave  effect.  Careful  tests  of  these 
ideas  can  be  made  by  using  the  data  obtzdned  in  the  present  work  via  linear  analyses 
and  direct  numerical  simulations. 
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APPENDIX  A 


RAPID  DISTORTION  THEORY 

Rapid  Distortion  Theory  (hereafter  RDT)  is  applied  to  study  the  response  of 
turbulence  quantities  to  one-dimensional  compression  due  to  a  shock  wave.  The 
state  of  homogeneous  turbulence  changes  significantly  when  it  is  subjected  to  mean 
strain,  but  the  nonlinearity  of  the  goverr  ng  equations  makes  it  impossible  to  de¬ 
velop  a  rigorous  theory  of  turbulence  under  straining  process.  For  cases  where 
nonlinear  effects  are  not  significant,  however,  one  can  solve  the  exact  linearized 
equations. 

When  the  time  scale  of  turbulence,  is  long  compared  to  that  of  the  mean 

deformation  (ISjq^/c  >>  1),  the  turbulence  has  no  time  to  interact  with  itself. 
Thus,  we  need  not  consider  the  nonlinear  terms  in  the  governing  equations  involving 
products  of  fluctuation  quantities,  and  we  can  obtain  the  linear  RDT  equations. 
The  viscous  terms  are  linear  and  can  be  included  in  the  analysis,  but  are  often 
neglected  and  will  be  here  for  simplicity.  RDT  for  one-dimensional  compression  is 
applied  to  study  the  response  of  turbulence  during  its  passage  through  a  shock  wave. 
In  the  following,  a  brief  review  of  the  RDT  procedure  is  given  for  one-dimensional 
compression  [Lee  and  Reynolds  1985,  Lee  1989]. 

We  consider  turbulence  subject  to  a  rapid  irrotational  strain,  where  density  is 
uniform  throughout  the  field  but  allowed  to  vary  in  time.  Since  the  mean  flow  is 
assumed  to  be  irrotational,  we  use  the  principal  axes  of  strain  rate  tensor  as  axes 
of  reference  so  that  the  mean  velocity  field  is 

f^adn  =  Sa0{t)  =  ra{t)S^0.  (.4.1) 

Using  the  equation  of  continuity,  we  can  then  express  tlie  evolution  of  density  p{t) 
in  the  form 


P(0  =  Pocxp 


(^4.2) 
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where  po  —  p(0)  and  the  trace  of  mean  strain  rate  tensor  5^,  signifies  the  expansion 
rate  (or  the  compression  rate)  of  the  mean  flow  field. 

The  turbulent  momentum  equation  in  the  RDT  is  given  by 


X  J1 

dt 


-SijUj 


1  dp 
pdx^' 


(/1.3) 


and  the  equation  of  continuity  is  unchanged.  The  dynamical  equation  for  the 
fluctuating  vorticity  ui  becomes 


X  TT 

dt  ^^^dXj 


(AA) 


where 

^ij  —  ■“  (^-5) 

is  the  deviatoric  component  of  the  mean  strain  rate  tensor.  The  first  term  on  the 
right  hand  side  of  (-4.4)  represents  the  production  of  vorticity  due  to  stretching  by 
the  incompressible  mean  strain  rate  S*j.  The  second  term  indicates  reduction  of 
vorticity  by  isotropic  expansion  (or  increase  by  compression). 

We  impose  periodic  boundary  conditions  and  represent  homogeneous  turbulence 
in  terms  of  Fourier  series.  The  dependence  of  Uj  on  the  position  x  in  (i4.3)  and 
(i4.4),  however,  poses  problems  with  this  approach.  To  remove  this,  we  use  a 
deforming  coordinate  system  [Rogallo  1981]  which  follows  the  mean  flow: 


f 

sa  —  , 

ea 


t, 


(4.6) 


where  the  factor  Cn  is  the  total  strain  in  the  o-direction  defined  as 


^oit) 


-  exp 


Jo 


In  the  deforming  coordinates,  the  vorticity  equation  becomes 


du>Q 

"dr 


(4.7) 


(4.8) 
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where 


To  —  Ta  —  To,  To  —  S{i  —  Fi  +  r2  +  r3 

are  the  reduced  strain  rate  and  the  dilatation  rate,  respectively.  The  solution  of 
(j4.8)  is 


P{t)  Po 


or  u>a(6'^)  =  eQC<;a(^,0). 


(A.IO) 


Here,  u;q(^,0)  is  the  initial  value  of  u;a(^,T),  and 

ea  =  — ,  =  €16263  =  — 

€0  P 

are  the  reduced  toted  strain  in  the  a-direction  and  the  dilatational  total  strain, 
respectively.  Note  that  Co  =  1  for  an  incompressible  mean  flow  field. 

Since  we  have  the  solution  (A. 10)  for  the  vorticity  evolution  in  an  explicit  form, 
it  is  easy  to  obtain  the  history  of  turbulence  statistics  in  terms  of  the  initial  values 
and  total  strains.  Thus,  the  evolution  of  the  vorticity  correlation  tensor  at  time  r 
is 


Vapi-r)  =  eaCpV^pW,  (A12) 

where 

V,p{0)  =  :j'cU,0)u;p{(,0).  (A.13) 

The  velocity  field  can  be  deduced  from  the  vorticity  field.  In  the  transformed 
coordinates,  the  Poisson  equation  for  the  velocity  is 


-ui  = 


1  d>jL;2 
€3  ^^3 


1  dui^ 
62  0^2 


(AM) 


where  the  transformed  Laplace  operator  is 


\  d-  ^  \  1 

eld(l  flcKj 


(.1.15) 
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The  equations  for  U2  and  U3  can  be  obtained  by  permuting  the  indices.  The 
solution  is  obtained  using  Fourier  expansions, 


«;n 

U,  =  ^  (A. 16) 

where  Ka  is  the  Fourier  wave  number  vector  in  the  deforming  reference  frame, 
'^Q  =  Ca^o  (^Q  is  the  wave  number  vector  in  the  undeformed  reference  frame). 

The  solution  for  uj  is 


ui 


V  ea  £2  / 


(yl.l7) 


where 


o 

X 


The  other  components  can  be  found  by  permutation  of  the  indices. 
Using  (.4.17),  we  find  the  velocity  spectrum  in  /c-space  to  be 


£:ii(/c,t) 


1 


(  —  )‘//22(«,t)  f  (--)^//33(K,r)  -  2(  — )( —  )//23(«, r) 

C3  ^2  e2  63 


(>1.18) 

where  H^j{k,t)  is  the  vorticity  spectrum  in  K-space,  and  the  relation  7723  —  ^32 
is  used.  The  vorticity  spectrum  at  a  later  time  is  obtained  in  terms  of  its  initial 
spectrum  as 


(>1.19) 
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The  initial  vorticity  can  be  expressed  in  terms  of  initial  velocities  in  the  Fourier 
space  as 

w/(/c,0)  =  -i€/^„/CmSn(K,0).  (A.20) 

Therefore,  we  can  represent  the  initial  vorticity  spectrum  Hij(K,0)  in  terms  of  the 
initial  velocity  spectrum  £jj(K,0),  for  example, 

/fll(/c,0)  =  /C2£33(/c,0)  +  /C3fJ22('<^>0)  —  2«:2'^3^'23('«>  0))  (>1.21) 

where  the  relation  JEJ23  =  £32  was  used. 

Using  the  relations  ( A.17),  (>1.19),  and  (A.20),  we  can  obtain  the  evolution  of 
the  velocity  spectrum  during  a  rapid  irrotational  strain. 

As  a  special  case  of  irrotational  strain,  we  consider  one-dimensional  compression 
in  the  ii-direction,  where 

Po  . 

€l  =  — ,  €2  =  63  =  1, 

P 

and 

€\  =  1,  ^2  =  €3  =  Cj  (A. 23) 

Vorticity  components  in  the  X2  and  X3  directions  are  amplified  by  a  factor  of 
density  increase  during  the  rapid  compression,  while  vorticity  in  the  xj-direction 
stays  the  same.  These  relations  hold  regardless  of  the  specific  form  of  the  initial 
spectrum.  In  general,  amplifications  of  Reynolds  stresses  do  depend  on  the  shape 
of  the  spectrum.  However,  for  an  isotropic  initial  turbulence,  amplifications  of 
Reynolds  stresses  are  again  independent  of  the  initial  spectrum:  for  example,  to 
compute  Uj  ,  triple  integration  of  (A. 18)  in  k  space,  using  (A.20)  and  a  general 
form  of  the  vorticity  spectrum  for  an  isotropic  turbulence, 

//,j(/c,0)  (A. 24) 

gi  VPS 


(A.22) 
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I  2 

Ui  = 


fOO  fOO  fOO 

/  /  /  Ell(K,T)dKldK.2dK2 

J  -oo  J  ~oo  J  ~oo 

1  —  k1)  +  /C9(l  —  /C?)  +  2/Co/Co 

=  —  /  E{K,0)dK  /  /  ^  sin  0d<pde, 

Jo  Jo  Jo  [(«l/ei)2+K2  +  «2j2 

(^.25) 


where 


/Cl 


/I  ~  '^2  .  o  ,  ~  «3  .  o  .  , 

—  =  cost/,  «2  =  —  =  sin  a  cos  (p,  ^^3  =  —  =  sine/ sin  0. 

K  K  K 


The  result  is  independent  of  /c.  Note  that  J  £'(/c,0)d/c  =  where  fJ(/c,0)  is 

the  initial  energy  spectrum  function  and  =  a'u'  at  <  =  0. 

The  same  result  can  be  obtained  by  properly  scaling  the  RDT  result  for  tur¬ 
bulence  under  incompressible  axisymmetric  expansion.  The  one-dimensional  com¬ 
pression  can  be  decomposed  into  isotropic  compression  ij/3  (or  density  change) 
and  incompressible  axisymmetric  expansion  S*^  as 


Sij  —  +  5* 


/ri/3 

0 

/2ri/3 

0 

0  \ 

ri/3 

0 

+  0 

-ri/3 

\  0 

0 

r,/3y 

\  0 

0 

-ri/3/ 

The  total  strain  can  be  accordingly  decomposed: 


e^ei'V  =(p^/p)l/3e*, 


where 


♦  /  *  *  ♦  \ 

e  ={ej,e.^,e3) 

is  the  incompressible  total  strain  vector  (i.e.,  ejC^Cj  =  1)  and 

fo  ~  ci«’2<^3  -  PoIp 


(A.26) 


(A.2T) 

(A.28) 


(A.29) 


144 


is  the  total  dilatation. 

Using  (i4.10)  and  (>1.20),  the  Fourier  amplitude  of  velocity  fluctuations  under 
one-dimensional  compression  (>1.17)  can  be  expressed  as 


1 

^3  2  , 

^2  2^ 

1  -  ^3 

€2 

9 

1  — «2  + 

—  «3 

Uio - KlK2U2o  - 

eoX 

L\e2 

^3  J 

^2 

^3  J 

(A.30) 


where  ''2/^2  initial  Fourier  amplitude.  Expres¬ 

sions  for  the  other  two  components  are  obteiined  by  cychc  interchange  of  indices. 

By  using  the  decomposition  (>1.26),  the  Fourier  amplitude  of  turbulent  velocity  in 
(>1.30)  can  be  expressed  as  a  product  of  the  isotropic  dilatation  and  incompressible 
contributions: 


u(/c,e)  =  {p/po)^^^u*{K,e*), 

where  from  (>1.30),  the  incompressible  part  u*  is  given  by, 


(A.31) 


q(«,c*) 


( "^4 +  ^^4)^10 - 

\  €2  / 


®3 

-T«1«2W2o 

^2 


^2 

-i^lK^u^o  , 

^3 


(>1.32) 


with  +  '^3/^3^*  The  appearance  of  the  total  dilatation 

as  a  common  factor  in  all  components  of  the  velocity  amplitude  implies  that  the 
anisotropy  in  a  turbulent  flow  is  not  affected  by  isotropic  compression  of  the  mean 
flow. 

The  energy  spectrum  tensor  after  one-dimensional  compression  is  expressed  as 

E,j(«,e)  =  {plpof^'^E*j{K,e*),  (A. 33) 

where  E*^  is  the  energy  spectrum  tensor  of  turbulence  subject  to  incompressible 
axisymmetric  compression. 

The  RDT  result  for  turbulence  under  one-dirnensional  compression  can  be  ob¬ 
tained  by  scaling  the  RDT  result  of  turbulence  under  incompressible  axisymmetric 
expansion  (see  (A.33)). 


14.5 
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APPENDIX  B 


LINEAR  INTERACTION  THEORY 

Some  aspects  of  the  interaction  of  turbulence  with  a  shock  wave  are  amenable  to 
linear  analysis.  A  fluctuating  quantity  in  a  compressible  flow  can  be  decomposed 
into  acoustic,  entropy,  and  vorticity  waves  [Kovasznay  1953].  Linear  interaction 
analysis  (hereafter  LIA)  of  a  shock  wave  with  various  linear  plcine  waves  were  per¬ 
formed  by  Ribner[1953],  Moore[1953],  Kerrebrock[1956],  Chang[1957],  and  McKen¬ 
zie  and  Westphal[l968].  In  general,  whenever  any  one  of  these  waves  passes  through 
a  shock,  it  generates  the  other  two  waves  downstream.  In  this  work,  we  follow  the 
methodology  of  Ribner,  where  the  main  interest  is  confined  to  the  interaction  of 
vorticity  waves  with  a  shock.  In  the  following,  a  brief  description  of  Ribner’s  anal¬ 
ysis  is  given  first,  followed  by  its  application  to  the  analysis  of  shock-turbulence 
interaction. 

B.l  Description  of  LIA —  Ribner’s  Analysis 

Ribner  formulated  the  interaction  of  a  plane  vorticity  wave  with  a  shock  wave 
as  a  boundary-value  problem  for  the  flow  in  the  region  downstream  of  the  shock. 
The  governing  partial  differential  equation  for  small-perturbation  rotational  flow 
was  derived  as  an  extension  of  Sears’  work  [1950],  boundary  conditions  on  the 
velocity  components  just  behind  the  shock  were  obtained  from  the  obhque-shock 
relations,  and  finally  the  rotation  term  in  the  governing  equation  was  evaluated  in 
terms  of  gradients  of  entropy  and  total  enthalpy  with  the  use  of  the  entropy  change 
across  the  shock.  The  initially  unknown  distortion  of  the  shock  wave  was  taken 
into  account  in  the  boundary  conditions  and  the  rotation  term  by  assuming  it  to 
be  sinusoidal  with  initially  undetermined  amplitude  and  phase. 

B.1.1  B’ormulation  of  the  Boundary- Value  Problem 

The  inclined  plane  sinusoidal  shear  wave  is  schematically  shown  in  F'igure  B.l. 
The  flow  is  viewed  in  a  plane  perpendicular  to  the  shock  and  to  the  wave  fronts. 
The  Weave  is  supposed  to  be  converted  downstream  by  the  mean  flow  with  velocity 
(most  variables  in  this  Appendix  are  used  exclusively,  and  those  commonly  used 
are  listed  in  the  Nomenclature  section).  The  pa.ssage  throiigh  tlie  shock  is  evidently 
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an  unsteady  process,  since  the  intercepts  of  the  inchned  shear  wave  with  the  shock 
wave  move  downward  along  the  shock  front. 

A  plane  oblique  sinusoidal  shear  wave  may  accompany  a  perturbation  velocity 
component  normal  to  the  plane  of  the  figure.  This  velocity  component,  however,  is 
parallel  everywhere  to  the  shock  and  is  unaffected  as  the  shear  wave  passes  through; 
this  component  also  has  no  effect  on  the  shock  wave  and,  therefore,  is  excluded  in 
the  following  analysis  without  any  loss  of  generality.  In  the  following  analysis,  the 
plcine  shear  wave  is  assumed  to  propagate  in  the  a;ia:2-pl3^ne. 

If  an  observer  moves  downward  along  the  shock  with  a  speed  V,  the  flow  has  an 
apparent  upward  velocity  component  V.  The  observer  speed  V  is  to  be  chosen  such 
that  the  resultant  mean  velocity  (in  the  observer’s  frame  of  reference)  is  aligned 
w'ith  the  velocity  in  the  disturbance  wave,  that  is,  V  =  tan  The  process 
appears  to  the  observer  as  an  interaction  of  a  steady  sinusoidal  shear  flow  with  an 
oblique  shock  wave.  Thus,  by  properly  choosing  the  frame  of  reference,  the  original 
unsteady  flow  problem  has  been  reduced  to  an  equivalent  steady  flow  problem. 

The  analysis  is  aimed  at  calculating  the  flow  field  downstream  of  the  shock  pro¬ 
duced  by  the  passage  of  a  sinusoidal  shear  flow  through  ihe  equivalent  oblique 
shock.  Because  of  the  nonuniform  upstream  velocity  field,  the  shock  wave  is  cor¬ 
rugated  and  introduces  vorticity  downstream  of  the  shock  wave.  If  the  upstream 
disturbance  wave  is  weak,  the  downstream  velocity  perturbations  are  also  weak 
compared  to  the  mean  velocity,  so  that  a  linearized  treatment  of  the  flow  field  is 
feasible. 

Ribner  derived  the  governing  partial  differential  equation  for  small-perturbation 
compressible  rotational  flow  with  gradients  in  entropy  and  stagnation  enthalpy.  In 
the  transformed  coordinates,  or  in  the  observer’s  frame  of  reference  shown  in  Figure 
B.2,  the  governing  equation  is  expressed  as 


and  t/’r/n 

we  have 


(1  -  A/2)V'cc  +  -n  (B.l) 

—  Using  Crocco  theorem  [Thompson  1984], 


n 


W  V  Br/  Or] )  ' 
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where  (  is  the  distance  in  the  streamwise  direction,  t]  is  the  distance  in  the  per¬ 
pendicular  direction,  W  is  the  streamwise  velocity,  AI  is  the  corresponding  Mach 
number,  H  is  the  stagnation  enthalpy,  s  is  the  entropy,  T  is  the  temperature,  is 
the  vorticity,  and  ip  is  a.  perturbation  stream  function  defined  as: 

dip/drj  =  -(1  -  M‘^)di'/dC,  =  {B.2) 

where  and  Wx)  are  the  perturbation  velocities  in  and  r)  direction,  respectively. 
The  final  flow  pattern  depends  crucially  on  whether  W  is  subsonic  or  supersonic, 
which  depends  on  the  Mach  number  corresponding  to  and  on  the  wave  inchna- 
tion,  6. 

The  boundary  conditions  just  downstream  of  the  shock  are  obtained  by  applying 
the  Rankine-Hugoniot  relations  across  the  perturbed  shock  wave.  By  geometry 
shown  in  Figure  B.3,  the  mean  velocity  components  normal  and  tangential  to  the 
undisturbed  shocks  are,  respectively, 

U =  W cos  T  =  sin  6. 

The  shear  wave  provides  a  perturbation  wj^  to  and  causes  indirectly  a  pertur¬ 
bation  (T^zo)  to  the  shock  wave  inclination  which  is  defined  as 

tano-  =  ^ — , 

0X2 

whose  magnitude  is  yet  to  be  determined.  The  effect  of  cr  is  equivalent  to  an 
increment  in  6.  The  associated  perturbations  to  U and  V  are 


dU^  —  ( IF4  4- u’^)  cos(^  +  (j)  -  cos  dV"  =  ( IF4  f  4 )  sin(^  +  cr)  —  IF^  sin 

Assuming  that  cr  and  w^/Wj^  are  small,  we  can  express  the  perturbations  as 

dU ^  ~  cos  ^  —  crlF^  sin  dV''  ~  sin  ^  4  crU'^^  cos  {B.3) 

The  change  in  shock-normal  velocity  across  the  sliock  wave  is  given  by  Rankine- 
Hugoniot  relations  as 
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£4  ^ 

Assuming  the  upstream  sound  speed,  =  constant,  the  change  in  the  downstream 
shock-normal  velocity  due  to  the  change  in  the  upstream  shock-normal  velocity  is 
given  by 


dU 

U 


dUA 

Ua 


(BA) 


where  m  =  U a/ U . 

On  the  downstream  side  of  the  shock,  the  perturbation  of  the  velocity  in  the  ^ 
direction,  W  ~  U  cos  6'  V  sin  6' ,  is 


=  {U  +  dU)  cos(e'  +  ^r)  +  (F  +  dV)sin{e'  +(t)-W,  (5.5a) 

and  the  perturbation  velocity  in  the  tj  direction  is 

^{U  +  dU)  sm{0'  +  a)  +  (V  +  dV)cos{e'  +  a),  (5.56) 

where, 

6'  =  tan~  ^  (m  tan  ^).  (5.5c) 

Dividing  both  sides  of  (5.5a)  and  (5.56)  by  U,  expanding  trigonometric  identi¬ 
ties  using  the  assumption  of  small  perturbations,  and  relating  dU/U  and  dV/U  to 
waI^^a  by  the  use  of  (5. .3)  and  (5.4),  one  gets 


U 

Wjj 

Tf 


- tan  f  1  -  2^^ - cos  O’  +  (  tan  O'  -f  mcr  )  sin 

^^A  /  V  7  +  1  /  X^^A 


\Va  m 


tan  (  1  —  2— - m  I  sin  O'  +  (  — —  tan  $'  +  mcr  ]  cos  O'  ~  a  sec  O' . 

)  \  7  +  1  /  X^Va  ) 


(5.6) 
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This  is  the  general  form  of  the  boundary  conditions  for  the  governing  equation 
(5.1)  downstream  of  the  shock  wave. 

In  the  present  analysis,  the  perturbation  is  a  sinusoidal  disturbance  velocity 
(associated  with  the  incident  shear  v»'ave)  parallel  to  WJ^, 

—  =  6  cos  kTJji 

Wa 

where  k  is  the  wave  number  and  is  perpendicular  to  Wj^.  The  corresponding 
argument  for  the  refracted  shear  wave  involves  rj  and  an  altered  wave  number  k' . 
The  wave  length  of  the  disturbance  should  match  along  the  shock  front,  that  is, 
^2  =  ^2’  ^2  =  kcos9  and  ki^  —  k'  cos  O',  and  along  the  shock  cos  0  =  t/^/x2 

and  cosO'  =  rj/xo-  Therefore, 

kT]A  =  k'r) 

and 

^=ecosk'ri  (5.7) 

Wa 

along  the  shock.  By  geometry  (Figure  B.l),  kcosO  =  k' cos  O'.  Since  the  dis¬ 
turbance  is  sinusoidal,  the  shock  inclination  a  is  expected  to  be  sinusoidal.  For 
generality,  a  phase  shift  is  allowed  for,  so  that  a  can  be  assumed  to  have  the  form 

a  =  cos  k'j^  +  bg  s\n  k'r]).  (5.8) 

Substitution  of  (5.8)  into  the  general  form  of  the  boundary  conditioris  (5.6),  yields 
the  final  form  of  the  boundary  conditions.  The  boundary  conditions  imposed  by 
the  shock  wave  on  the  perturbation  velocity  components  parallel  to  (  and  rj,  nn- 
mediately  behind  the  shock  wave  are 


eU 


05  /  1  ^  A  ■  a'  ( ^ 

—  1  -  2 - m  “f  rn  sinW  -  1  -  2  — 

mV  7+1  y  V  7+1 


sin‘  O' 


)  l  »1U  f 

cosu  4  - .  COS  k  n 


COS  I 


^5  ^7  1 


■-  (  1  -  2- - rn  t  ]  sin  O' 

mV  7  I  1 


sin  k'r] 


respectively.  In  (Z?.8),  the  parameters  ag  and  65  governing  the  shock  inclination  (t 
are  undetermined. 

For  the  solution  of  the  governing  equation  (5.1),  the  vorticity  term  on  the 
right-hand  side  must  be  evaluated  for  the  region  behind  the  shock  wave.  Down¬ 
stream  of  the  shock,  the  stagnation  enthalpy  H  and  the  entropy  s  (and  hence  the 
vorticity)  are  constant  along  streamlines,  and  in  the  linear  theory  the  streamlines 
are  approximated  by  lines  of  constant  7/.  Thuj,  dH /dr]  and  ds/dr]  may  be  evaluated 
at  the  shock  and  the  result  also  holds  downstream  for  the  same  77. 

The  total  enthalpy  upstream  and  at  the  shock  is 


The  entropy  upstream  of  the  shock  is  constant  by  virtue  of  the  assumption  of 
constant  pressure  and  density  there.  The  entropy  change  across  the  shock  is  given 
in  terms  of  the  upstream  velocity  by 


[Thompson  1984].  DifTerentiating  and  expanding  the  above  expression  under  the 
assumption  of  small  perturbations  gives 

along  the  shock. 

Substitution  of  (5.10)  and  (fl.ll)  into  the  governing  equation  (5.1)  yields 


(1  -  4-  =  Urn 


2  COS  ^  d 


(r.) 


cos^  0  drj 

U  cos  $'{Tn  —  1)^-^  I 
or}  \ 


W. 


—  a  tan  6 


(5.12) 


where  the  right-hand  side  is  evaluated  along  the  shock  (rj  =  0)  and  expressed  as 
function  of  r]  alone. 

Substitution  of  (5.7)  and  (5.8)  for  Wj^/W a  and  a  for  the  sinusoidal  form  of  the 
upstream  disturbances  converts  (5.12)  into 


(1  -  M^)0^^  4  V>T,T, 
-  Uel-k' 


(  1 

sec  o'  4  2(m  —  1 )  cos  O'  4  05 - sin  O' 


m 


■  ,1  ,  ,'L  ~  ■  ni  ,1 

sin  k  Tj  4  k  - -  sin  a  cos  k  77 


m 


(5.13) 

where  the  relation  tan^^  =  Tnta.n0  is  used  to  eliminate  0.  Equation  (5.13)  is  the 
partial  differential  equation  governing  the  flow  downstream  of  the  shock  subject  to 
the  boundary  conditions  (5.9),  which  are 


dxi> 
Ot]  ’ 


Wr,  =  --(1 


M-) 


50 
5C  ■ 


1.53 


B.1.2  Solutions  to  the  Problem 

The  nature  of  the  governing  equation  (B.13)  depends  on  the  sign  of  (1  —  A/^): 
for  flows  with  M  <  \  the  partial  differential  equation  is  elliptic,  and  for  flows  with 
M  >  1  it  is  hyperbolic.  In  the  following,  therefore,  two  different  sets  of  solutions 
to  (B.13)  are  presented. 

Solution  when  Flow  Downstream  of  the  Shock  is  Subsonic,  M  <  1 

The  governing  differential  equation  (fl.l3)  can  be  written  in  abbreviated  form 
as 


—  -h'Ue{As  sink'r]  —  Bs  cos  k'r}),  (-5-14) 

where 


(rn  —  1 

Ag  =  sec  +  2(m  —  1 )  cos  O'  +  ag - sin  O' 


m 


m 

2 


(B.15) 


A  particular  solution  to  (5.14)  is 


/  If  ( ■  l'  u'  \ 

i'p  =  t/  €  I  —  sin  A:  7/  -  —  cos  A:  77  1 


(5.16) 


A  homogeneous  solution  is 


■/’(7  =  U  te 


((^  cos  S'— T7sin^')''os  S' 


>:  <  c  cos 


^  k'  cos  o' 


((  sin  o'  +  /3^77  cos  O' ) 


f  d'  sin 


k'  cos  O' 

”52 


-{(  s'm  o'  +  cos  O' ) 


(5.17) 


l.')4 


where  c'  and  <i*  are  the  constants  of  integration,  and  (3^  =  I  —  /a^. 

The  complete  solution  for  the  perturbation  stream  function  V*  =  V’p  +  i’C 
tains  four  undetermined  parameters,  as,l>s,c',  and  d!  which  are  to  be  determined 
by  matching  velocity  components,  =  dxl}/dr}  and  Wtj  —  —(i^dtp/dC,  at  the  shock 
wave  ((^  =  0)  using  boundary  conditions  (B.9). 


+  DsFs 

-  DsEs 

C|  +  l,|  ' 

”  Cl  +  Dl  ’ 

d'  - 

1  r\ 

k' (3^  cos  O' m 

(B.18) 

U  — 

where 


Cs  -  -  1)^  + 


2  ,  2(m  -  1) 

7  +  1 


sin  o'  cos  o' 


=  1  +  (m  -  l)cos^^'l 


r7  (  7  ~  1  \  ^  o' 

Es  =  2(^1-  +  2(m  -  1) - - 


Fg  =  [2(m  —  1 )  sin  cos  . 


(B.19) 


Solution  when  Flow  Downstream  of  the  Shock  is  Supersonic,  M  >  1 
The  particular  solution  to  (B.13)  is  the  same  as  that  for  M  <  \,  which  is 

i’P  =  sin  k  Tj, 

where  the  fact  that  the  final  solution  yields  65  =  0  (and  hence  B5  0)  is  used  to 
delete  Bg  term  at  the  outset. 

The  complementary  solution  satisfying  (B.13)  has  the  general  form  of 
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V-C  =  /(C  +  0M>r))  +  5(C  -  0w^), 

where  =  A/^  — 1.  The  function  /  represent  Mach  waves  inclined  downward  by  the 
Mach  angle  from  the  ^-axis  and  the  function  g  represents  Mach  waves  inclined 
upward  by  the  same  angle.  For  upstream  shear  waves  with  inclination  0  <  ^  <  7r/2, 
the  g-family  of  Mach  waves  represent  disturbances  overtaking  the  shock  wave  from 
behind  (propagating  upstream).  However,  its  propagation  toward  upstream  beyond 
the  shock  wave  is  unphysical,  since  the  Mach  waves  can  only  propagate  with  the 
speed  of  sound  of  the  medium  relative  to  the  mean  flow  speed  while  the  upstream 
mean  flow  speed  is  <dways  supersonic,  >  1:  Therefore,  ^-function  must  be 

zero  for  upstream  shear  waves  with  inclination  0  <  ^  <  7r/2.  For  the  same  reason, 
/-function  must  be  zero  for  upstream  shear  waves  with  inclination  k/2  <  6  <  ir. 
In  the  following,  the  discussion  is  limited  to  the  specified  range  0  <  ^  <  n  12  due  to 
the  symmetry. 

From  {B.9)  with  65  =  0,  the  function  /  reduces  to 

/  ~  sin{a(C  +  0wn)}- 

Along  the  shock  front,  where  (  =  TflanO', 

a(C  +  BwV)  = 

A  suitable  complementary  function  is  therefore 

Utc"  .  k'{(  + /SwV) 
wr  —  — rr-  sin  — - — , 

k'  +  tan  B' 

where  c"  is  a  constant  of  integration. 

The  complete  solution  for  the  perturbation  stream  function  is 


t/>  =  i/’/>  -f  V’c  ^ 


A5  sin  k'r]  -t-  c"  sin 


k'iC  +  0w^) 

i3yj  -f  tan  6' 


l.'>6 


(5.20) 


where  the  two  arbitrary  parameters  05  (occurring  in  Ag  of  (5.15))  and  c"  remain 
to  be  determined.  These  two  parameters  are  determined  by  matching  velocity 
components  derived  from  (5.20)  with  the  boundary  conditions  given  in  (5.9).  The 
resulting  expressions  for  the  parameters  are 


o-S 


C's  +  GsF's 


7^  (  m  ’ 


where 


(5.21) 


C'c  =  2- - m  —  2[1  +  (m  —  1 )  cos^  6'] 

^7  +  1 

D'^  —  (m  -  1)[1  -I-  (m  —  1)  cos^  $'] 

E's  =  {m  -  1  )'^  sin  $'  cos  tan  6' 

Eg  —  2(m  —  1 )  sin  6'  cos  $' 

^  +  tan^'  ""  -  ^')  (5.22) 

with  the  Mach  angle  =  cot“^5u>- 
B.1.3  Summary  and  Discussion 

The  interaction  of  a  plane  shear  wave  and  a  shock  wave  is  schematically  shown  in 
Figure  2.1.  The  method  of  predicting  the  flow  field  downstream  of  the  shock  wave, 
produced  by  convection  of  an  oblique  sinusoidal  shear  wave  through  the  shock,  has 
been  briefly  described.  In  this  section,  the  main  results  arc  summarized  in  more 
compact  form,  and  are  simplified  to  help  the  geometrical  interpretation. 

Criterion  on  M 

Although  the  streamwise  velocity  C  dowmstream  of  the  sjjecified  normal  shock 
is  always  subsonic,  the  nature  of  the  flow  depends  primarily  on  the  streamwise 
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velocity  in  the  transformed  frame  of  reference  (Figures  B.2  and  2.1),  which  may 
be  subsonic  or  supersonic.  Two  forms  of  the  solution  for  all  flow  quantities  thus 
appear:  one  for  the  subsonic  range  M  <  1  and  the  other  for  the  supersonic  range 
M  >  1.  Since  M  depends  on  the  initial  Mach  number  V aI^^A  the  inclination 

angle  6,  the  equation  for  the  dividing  line  A/  =  1  gives  a  relation  betw'een  the 
critical  value  of  6  and  the  upstream  Mach  number  UaI^a 


Oct  =  ±  tan' 


f{^  4-  l)(Tn  -  1) 

2m^ 


(i?.23) 


where  m  =  UaIU  is  determined  by  UaIo-a- 

For  waves  incident  at  angles  less  than  the  critical  angle,  a  pressure  wave  generated 
at  time  t  catches  up  with  that  generated  at  time  t  +  because  the  upstream 
shear  wave  travels  along  the  shock  wave  slower  than  the  pressure  wave  generated: 

tan^  0  <  a?  -  .  Note  that  the  upstream  shear  wave  travels  with  the  velocity 

U A  tan  6  along  the  shock  wave,  and  the  generated  pressure  wave  propagates  radially 
outward  with  the  speed  of  sound  a  relative  to  the  mean  flow  speed,  U.  In  this  case, 
pressure  waves  are  superposed  onto  each  other,  and  disable  each  pressure  wave 
from  propagating  independently,  resulting  in  an  exponential  attenuation  in  the 
shock-normal  direction.  For  waves  incident  at  angles  larger  than  the  critical  angle, 
the  pressure  waves  generated  at  different  times  propagate  independently  from  the 
others,  because  tan^  0  >  a?  —  . 

Downstream  Velocity  Field 
By  u.sing  (B.2),  for  M  <  1, 


— —  5  cos(A:(,(x2  ^  •^'1  tan  O')  A  >^3]  +  n(xi )  cos[^-2(x'>  —  tan  O")  +  ^p] 


; — -  /luTIfxi )  sin[A:2(3-2  '  tan^")  +  Sp\,  (B.24) 


where  \u'a\  -  's  the  amplitude  of  the  upstream  sinusoidal  velocity  and 
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n(ii) 


cost/V-5  >  -5 
m  /3 


=  tan 


=  —  tan  ^ 


t/^  tan  Q' 

/32 


,  -I  CS0W  -  dsiSinO' 

tan  — — - - 

dsPw  +  C5  tan  6' 


with  C5  =  c'k'0n,cos6'  ,ds  =  d' k' 0tu  cos  .  The  functions  >15  and  Bg  are  given  by 
(B.IS),  and  ag,bg,c',  and  d’  are  the  constants  of  integration  given  in  (5.18). 

For  M  >  1, 


=  5  cos  ^2(^2  “  ^1  tan  0')  +  n  cos  ^2(^2  ~  tan  0") 


- — —  =  /3u,n  cos  ko{x2  —  X]  tan  0"), 


(5.25) 


where, 


m 


^2  —  k'  cos  0'  :=  k  cos  0 

cosB  sinfixf 

rn  cos{0'  ~  1.1  \f  ) 


II 


0"  =  e‘  -  fiM 

fiM  =  cot"^ 


with  cg  =  c" j3u)  cos  0' .  The  function  A5  is  still  given  by  (B.15),  and  ag  and  c"  are 
evaluated  in  (5.21). 

The  cosine  in  the  first  terms  of  (5.24)  and  (5.25)  are  constant  adong  lines  X2  — 
xj  tan  =  constant,  which  are  inchned  at  an  angle  6'  with  the  horizontal  and  are 
thus  parallel  to  the  ^-axis.  Since  is  parallel  to  ^  and  Wjj  is  parallel  to  tj,  the  first 
terms  represent  contribution  to  the  vorticity  wave. 

The  remaiining  terms  in  (5.24)  and  (5.25)  involving  the  factor  IT  correspond  to 
an  irrotational  velocity  field,  or  potential  flow.  If  derivation  is  traced  backward, 
the  Il-terms  are  found  to  originate  from  the  complementary  solution,  which  is  a 
solution  with  zero  vorticity,  Q  =  0.  These  remaining  terms  in  (5.24)  and  (5.25) 
represent  contributions  from  pressure  waves.  For  the  case  M  >  1,  this  pressure 
wave  propagates  in  the  form  of  Mach  waves  with  the  speed  of  sound  relative  to 
the  stream  velocity.  For  the  case  A/  <  1,  the  resultant  pressure  pattern  does  not 
propagate  with  the  speed  of  sound,  but  it  can  be  represented  as  a  superposition 
of  cylindrical  sound  waves  which  individualb'  propagate  with  sonic  speed.  The 
resultant  perturbation  velocity  attenuates  exponentially  with  ihe  dist,.nce  from  the 
shock  wave,  xj  (see  (5.24)). 


Shock  Wave  Perturbation 

From  (5.8),  the  local  perturbation  in  the  shock  inclination  angle  can  be  repre¬ 
sented  as 

cr  =  e(ag  cos  k2X2  -I-  ^5  sm  ^’22^2 )» 

where  ag  and  bg  are  evaluated  in  (5.18)  for  Af  <  1  and  (5.21)  for  Af  >  1  {bg  =  0 
for  Af  >  1 ). 

At  any  point  along  the  shock,  X2,  the  local  shock  deflection  (  from  the  plane 
X|  =  0  is  obtained  by  integration  of  the  slope  <t  ~  d(,ldx2'- 
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^  =  /  (Tdx2 


^  Wa - '*■  (5.26) 

where  ^gj,  =  tan” ^(05/65)  is  the  phase  shift. 

For  a  given  wave  length,  the  amplitude  of  this  sinusoidal  corrugation  in  the  shock 
wave  is  proportional  to  the  factor  <^5  +  ^5  s^nd  edso  to  the  upstream  turbulence 
intensity. 

B.2  Application  of  LIA  to  Shock-Turbulence  Interaction 

Homogeneous  turbulence  can  be  represented  as  a  spectrum  of  waves  with  random 
orientations  and  wave  lengths.  In  a  solenoidal  (or  incompressible)  velocity  field,  a 
velocity  vector  in  Fourier  space  is  perpendicular  to  its  corresponding  wave  number 
vector  due  to  the  continuity  constraint.  (Any  incompressible  velocity  field  can  be 
decomposed  only  into  vorticity  waves.)  The  velocity  vector  corresponding  to  a 
vorticity  wave  in  the  Fourier  space  meets  this  constraint,  as  is  schematically  shown 
in  Figure  2.1.  When  turbulence  is  convected  into  a  shock,  the  individual  vorticity 
waves  are  abruptly  altered  producing  acoustic  and  entropy  waves  downstream  of 
the  shock.  This  process  is  schematically  sh  wn  in  Figure  2.1.  The  amplitude  and 
wave  length  of  downstream  waves  can  be  related  to  those  of  the  upstream  vorticity 
wave  through  linear  theory.  Therefore,  given  the  upstream  statistics  of  a  solenoidal 
turbulence  field,  we  can  obtain  its  downstream  turbulence  statistics. 

B.2.1  Turbulence  Modification 

An  upstream  vorticity  wave  is  associated  with  a  velocity  field  which  can  be 
decomposed  into  three  velocity  components  in  cylindrical  coordinates  as  shown 
in  Figure  2.1.  Each  velocity  component  can  be  represented  in  terms  of  Fourier 
coefficients  as 


diif  =;  du^  exp(ik  •  x), 


(5.27) 
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wlierc  the  wave  number  vector  k  lies  in  the  a-i,r-plane,  making  an  angle  9  witli  the 
r-axis.  'The  velocity  t)f  the  refracted  vorf icily  wave  can  also  be  t'xpressed  as 


(lu[  (/i7j  exp(tk^  •  x),  (/^.28) 

where  k*  is  the  new  wave  number  vector,  making  an  angle  9'  wdth  the  r-axis.  The 
radial  components  of  k*  and  k  are  etpial,  and  the  dependence  of  k^  on  k  is  expressed 
through  the  dependence  of  9'  on  9: 

9'  tan^  *(;/(  tan  0),  {B.29) 

where  -  J  *  >  1  is  the  density  ratio  across  the  shock  wave. 

Amplitudes  of  downstream  velocity  components  associated  with  the  vorticity 
wave  can  l)e  related  to  those  of  ui)stream  velc)city  comj)t>nent.s  as 


where  k”  is  the  wave  number  vector,  making  an  angle  0"  with  the  r-axis.  Again 
the  radial  component  matches  that  of  k.  The  amplitude  of  downstream  velocity 
components  associated  with  the  acoustic  wave  can  be  related  to  those  of  upstream 
velocity  components  as 


where 


du^  =  —  du], 
du^  =  r  dur, 

d^;^  =  0,  (5.33) 


_  _  cos  0'  -  sin  0' 

Cos  0 


r  = 


sin  0'  +  cos  0' 

sin  0 


(5.34) 


with  n  =  1  for  0  <  d  <  ^cr  and  n  =  0  ior  0ci  ^  0  ^  1^12.  That  is,  the  transfer 
functions  in  (5.31)  and  (5.34)  are  valid  for  both  subsonic  and  supersonic  cases 
with  the  appropriate  choices  of  8g  and  8p. 

Turbulence  components  in  the  Cartesian  coordinates  are  obtained  using  the  re¬ 
lations 


duo  =  dur  cos(f)  —  du^  sind>,  du^  =  dur  s\u<f)  -j  da^  cos4>,  (5.35) 

and  conversely 

duj.  ■—  duocostf)  +  du2s'u\<p,  du^  —duoslmi)  du^cosd.  (5.36) 

The  same  relations  apply  to  the  downstream  velocity  comj)onents  du'  and  du" . 
Using  (5.30)  and  (5.33)  -  (5.36),  we  can  represent  downstream  waves  in  terms  of 
the  amplitude  and  wave  number  of  the  corresponding  upstream  wave.  Upstream 
and  downstream  vcirticity  waves  represented  in  (5.27'  and  (5.28)  can  be  related 
as 
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(5.37) 


du^  =  Xdui, 

du2  =  Ydur  cos<p  —  du^siTK^^ 
diL^  —  ^  dur  sm4>  -f  du^  cos(f>. 

The  velocity  fields  of  the  corresponding  acoustic  wave  downstream  of  the  shock  is 
expressed  in  (5.32),  w’here 

Hduj , 

Tdur  cos4>, 

Tdur  sin<f).  (5.38) 

The  quantity  u^{x)  can  be  obtained  from  the  Fouriei  transform  as 


j-'" 
auj  = 

du2  = 


u 


x) 

•  x)dui(k), 


(5.39) 


where  the  triple  integration  is  extended  over  (  —  00,00)  in  each  component  of  k  = 
(ki,k2,k2).  By  defining 


we  can  rewrite  (5.39)  as 


dii,(k)  =  uj(k)dk, 


u,(x) 


J  u,(k)  exp(ik  •  x)dk. 


(5.40) 


(5.41) 


Tf  th#'  fi'-ld  is  homogeneous,  the  two-point  correlation,  u,(x;)iij(x')  can  be  repre¬ 
sented  as 


u,(x)uj(x') 


II 


exp[i(k-x  -f  Tx^)]  u,(k)uj(lj  dkdl, 


(5.42 
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where  the  triple  integrations  are  extended  over  (  —  00,00)  in  each  component  of  k 
and  1,  and  the  overhne  (•)  means  the  average  over  an  ensemble.  By  the  use  of  the 
orthogonality  of  the  Fourier  modes,  {BA2)  can  be  rewritten  as 


ui{x)uj{x')  =  y  exp(ik  •  r)  ui(k)u*(k)dk,  (B.43) 

where  r  =  x'  —  x,  u*(k)  is  the  complex  conjugate  of  Uj(k).  The  energy  spectrum 
tensor  JEJjj(k),  the  spectral  density  of  UjUj,  is  defined  as 

£;.j(k)  =  £.(k)u*(k)  (B.44) 

so  that 

ujuj  =  J  Eij{k)dk.  (5.45) 

In  order  to  predict  the  changes  in  the  second-order  turbulence  statistics  through 
the  interaction,  we  multiply  both  sides  of  (5.37)  by  their  complex  conjugates  and 
add  the  last  two,  resulting  in 

du^du^*  =  |X]2  du^dul 

du2du2  +  du^du^*  =  \Y\^  durdu*  du^da*^.  (5.46) 

From  geometry  (see  Figure  2.1), 


dur  —  dui  tan^, 

du*  —  du^  tan^,  (5.47) 

and  also  from  the  coordinate  transformation  (5.35)  we  have, 

durdu*  -f  dUf^du*^  =  du-ydu^  +  du^du^.  (5. 48) 

Thus,  using  relations  (  5.40),  ( 5.45),  and  (  5.48),  the  enseml)le  average  of  (  5.46) 
can  be  expressed  as 
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u\u\^dk'  =  \X\r  iliSitfk, 


{u'^u'2*  +  u'2u'^*)dK  =  (|y|  —  1)  tan'^^  uiUjffk  +  {u2^2  U2U^)dk.  (J5.49) 


If  operations  similar  to  those  performed  on  {B.37)  are  applied  to  (5.38),  we  get 


u"u'^*dk"  =  |H|^  iilSjdk, 

{U2U2*  +  u'^u'^*  )dk"  =  irf'^  tan^^  uj  Ujdk.  (5.50) 

The  mean-square  velo'^ity  components  associated  with  vorticity  waves  follow 
directly  from  integration  of  the  spectral  density.  Integration  of  both  sides  of  (5.49) 
yields 


.Yl'^Eiilk)  dk, 


ii  2  ^  3  —  u  j  -f‘ 


ul  +  J 


(ir|'  -  l)tan2^  5n(k)  dk. 


(5.51) 


Similarly,  integration  of  (5.50)  yields  the  mean-square  fluctuating  velocity  compo¬ 
nents  in  acoustic  waves: 


=  5n(k)  dk, 

n^tan^d  En(k)  dk.  (5.52) 

Thus,  the  mean-square  velocity  fluctuations  behind  the  shock  are  given  in  terms 
of  those  ahead  of  the  shock,  the  single- wave  transfer  functions  A'’,T,  5,  and  T, 
and  the  longitudinal  spectral  density,  Eii(k).  Note  that  the  single-wave  transfer 
functions  are  functions  of  k  and  0,  where  k  |k|. 
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Due  to  the  interaction  of  turbulence  with  a  shock  wave,  aerodynamic  noise  is 
produced  behind  the  shock  wave  in  the  form  of  fluctuating  pressure,  p’ .  The  fluc¬ 
tuating  acoustic  pressure  is  related  to  the  velocity  fluctuation  as 


or 


p  =  —pWwp 


.D 


-7A/ 


2^ 

W' 


(5.53) 


where  wp  is  the  (^-component  of  the  perturbation  velocity  associated  with  the 
pressure  fluctuation,  and  p^  is  the  mean  downstream  pressure.  (5.53)  can  be  rec¬ 
ognized  as  the  linearized  Bernoulli  equation  applied  to  the  velocity  in  the  acoustic 
wave. 

Substituting  for  M  and  IT,  and  using  (5.24)  and  (5.25)  for  wp  without  vorticity 
wave  contribution,  (5.53)  results  in 


-,D 


\'wj[\  27m  n 

(7  4-  1  )rn  —  (7  —  1 )  cos  O' 


cos[A:2(x2  —  xi  tan  0^^) -f- <5p],  (5.54) 


where  downstream  pressure  wave  is  expressed  in  terms  of  the  shock  strength  for 
the  corresponding  upstream  vorticity  wave  (17^  =  U\).  (5.54)  can  be  interpreted 
as  the  relation  between  the  Fourier  coefficients  of  the  downstream  pressure  wave 
and  upstream  velocity  fluctuations: 


ui/  cos  0  2jm  n 

~  (7  +  l)m  -  (7^T] 


(5.55) 


where  p  and  il]  are  the  Fourier  coefficients  of  the  pressure  and  streainwise  velocity 
fluctuations,  respectively.  Pressure  fltictuations  can  be  calculated  using  the  relation: 


/ 


p(k)^(Y)d^k, 


(5.56) 


where  ^  is  the  complex  conjugate  o  fp. 
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Following  a  similar  procedure,  we  can  relate  the  downstream  power  spectra  to 
the  upstream  power  spectra,  and  all  their  derived  quantities  of  interest  (such  as  two- 
point  correlations,  turbulence  length  scales,  and  dissipation  rate)  to  their  upstream 
counterparts. 

B.2.2  Shock  Front  Distortion 

The  local  perturbation  in  the  shock  front  inclination  angle  a  in  {B.8)  can  be 
considered  as  the  shock  inclination  when  r-axis  in  Figure  2.1  coincides  with  X2- 
axis.  In  general,  the  shock  wave  inclinations  in  xiT2  and  rir3-planes,  a2  and  <73 
respectively,  can  be  expressed  as 

a2  -  cr  cos  4>  a.nd  <73  =  cr  sin  4>. 

By  the  use  of  (5.8),  a  —  (f72,cr3)  can  be  written  as 


«*xi)[i(kft -x/i  +  (5.57) 


w 


here 


cr' 


(5.58) 


-=  {d/dx2,d/dx2),  k/j  “  (^'2i^'.3)>  ^  is  local  shock  front 

displacement  from  the  mean  shock  position;  f^i(=  U is  the  mean  upstream  flow 
speed. 

For  isotropic  turbulence,  v.iric.nces  of  the  local  inclination  angles  are  related  by: 


■'  •>  2 

rr~  <72  +  <73 

2<t|.  (5.59) 

Using  (  5..57 ),  ( 5..58 )  and  (5.59),  vve  can  express  the  variance  rr.;  for  an  isotropic 
upstream  turbulence  as 


(5.60) 


where  ujUj  =  5ii(k)  and  triple  integration  is  extended  over  the  wave  number 
space.  We  can  evaluate  this  integral  numerically  in  the  spherical  coordinate  system 
defined  as 


ki  =  A:  sin^,  k2  =  kcosO  cos4>,  k^  =  k  cosO  sin<f>. 


and 


dk  =  k^  cosd  d<f)d9dk. 

In  this  coordinate  system,  5ii(k)  can  be  expressed  as 


Eik)  o 

Eu{k,d,4>)  =  cos  e, 


where 


/ 


T  3  2 
E{k)dk  -ul, 


and  Uo  is  the  rms  of  the  velocity  fluctuations  in  one  direction. 
Now  (B.60)  can  be  rewritten  as 


(5.61) 

(5.62) 

(5.63) 


(5.64) 


- -  [  f  f  — cos*"^  (oc  4  be)  k"  cos9d4>d9dk 

2U\  Jo  Jo  Jo  47rA:2  ' 

— ^  /  E{k)dk  f  (o^  4  6^)  cos^9d9 
2Vt  Jo  Jo 


(4  -  4) 


(5.65) 
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where  symmetry  in  6  is  used  in  the  eveiluation  of  the  integral.  Note  that  the 
result  of  the  integration  is  independent  of  the  specific  form  of  the  upstream  energy 
spectrum,  E{k).  Since  the  remaining  integral  is  only  a  function  of  the  upstream 
Mach  number,  the  variance  of  the  shock  inclination  angle  is  dependent  only  on 
turbulence  intensity  and  the  Mach  number  upstream  of  the  shock  wave. 

Using  (B.57),  the  local  shock  displacement  ^  in  (B.26)  can  be  expressed  as: 

{  =  t^|5|  exp[i{kh  ■  Xh  +  ^sh)l-  (BM) 

ikh 

The  local  shock  front  curvature  in  the  12  direction,  K2  —  d^^/dx^,  can  be  ob¬ 
tained  by  differentiation  of  (5.57)  as 

=  exp[i(k/, -x/j  +  (5.67) 

kh 

Using  (5.58)  and  (5.66),  and  following  the  same  procedure  which  was  used 
to  determine  the  variance  of  the  local  shock  inchnation  angle,  we  can  obtain  the 
variance  of  the  local  shock  displacement: 


_  1  too  /"JT  y*2)r  fj/U)  1 

^  II  ~To  +  ^9)  k^  cos0d4>d0dk 

U{  Jo  Jo  Jo  47rA:2  ^ 

-  r  r'\ai + 4)  cosws. 

U,  Jo  Jo 


(5.68) 


Using  (5.58)  and  (5.67),  .we  find  the  variance  of  the  shock  front  curvature  in 
the  X2  direction  to  be 


1  E(k) 

- 9  III  — 79  cosJ0  ( A.'"  cos^^ cos‘^0)  (a|-  +  6^)  k^  cos0dd>d0dk 

2f7j“  Jo  Jo  Jo  47rfc‘- 


.3 

VeF^ 


f  k'^Eikjdk  f 

Jo  Jo 


r/2 


/  2 
(«v  ^ 


h"^ )  cos^0d0. 


(5.69) 
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Note  that  variances  of  shock  front  displacement  and  shock  front  curvature  are 
dependent  on  the  shape  of  the  upstream  turbulence  spectrum  as  well  as  on  the 
mean  flow  Mach  number  and  turbulence  intensity  upstream  of  the  shock  wave. 
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Shock 


I*  IGURE  B.  1 .  Convection  of  plane  oblique  sinusoidal  shear  wave  through  the  shock 
wave. 
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Shock 


Figure  B.2.  Symbols  and  coordinate  systems  used  in  LI  A 


Figure  B.3.  Geometrical  relations  across  the  shock,  with  and  without  perturba¬ 
tion  <T  in  the  shock  inclination  angle. 
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APPENDIX  C 


REMOVAL  OF  ALIASING  ERRORS 
IN  SIMULATION  OF  COMPRESSIBLE  FLOWS 


In  computational  fluid  mechanics  aliasing  errors  arise  from  the  non-linear  terms. 
The  product  of  two  quantities  resolved  on  a  mesh  of  size  N  results  in  a  function 
with  higher  frequency  components  that  can  not  be  resolved  on  the  mesh.  These 
components  are  then  “aliased”  back  in  the  computational  mesh  and  contaminate 
the  solution.  Aliasing  errors  have  been  reported  to  lead  to  inaccurate  simulations 
[Kim  et  al.  1987,  Spalart  1988]. 

In  a  simulation  of  incompressible  Navier-Stokes  equations,  aliasing  errors  are  pro¬ 
duced  from  the  product  of  two  velocity  components  in  the  convective  term.  (For 
simplicity,  we  only  consider  the  one-dimensional  problem  and  choose  the  Fourier 
collocation  method  as  the  numerical  differentiation  scheme.)  A  real  function  de¬ 
fined  on  N  grid  points  can  be  represented  with  N/2  complex  Fourier  coefflcients. 
The  multiplication  of  two  variables  generates  contributions  up  to  the  wave  num¬ 
ber  N  which  is  beyond  the  resolution  limit.  These  contributions  lead  to  aliasing 
errors.  Since  we  know  the  exact  origin  and  destination  of  aliasing  errors  we  can 
remove  them.  First,  we  expand  the  wave  number  domain  from  {  —  N/2  -|-  I,  N/2) 
to  {  —  3N/A  -f  l,3V/4)  with  zero  Fourier  coefficients  for  the  expanded  modes.  Mul¬ 
tiplication  is  then  performed  in  the  physiccil  space  with  the  expanded  grid.  The 
result  is  then  Fourier  transformed  and  the  additional  modes  are  eliminated. 

In  a  simulation  of  compressible  Navier-Stokes  equations  using  conservative  vari¬ 
ables,  aliasing  errors  arise  from  both  multiplication  and  division  operations.  Divi¬ 
sion  is  necessary  to  get  the  velocity  components  and  temperature  from  the  conserva¬ 
tive  variables.  Removal  of  aliasing  errors  stemming  from  multiplication  is  possible 
in  an  analogous  manner  as  in  the  incompressible  case.  However,  w’e  can  not  remove 
all  aliasing  errors  explicitly,  since  division  distributes  the  aliasing  error  throughout 
the  wave  number  domain. 

Good  spatial  resolution  is  helpful  to  minimize  the  effect  of  abasing,  since  the 
aliasing  errors  are  not  significant  for  a  well-resolved  flow  [Canuto  et  al.  1988]. 
Using  a  numerical  scheme  which  conserves  physically  important  (jnantities,  such 
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as  mass,  kinetic  energy,  and  total  energy,  helps  to  control  the  instability  resulting 
from  aliasing  error  [Feiereisen  et  al.  1981,  Blaisdell  et  al.  1990]. 

In  the  following,  we  propose  a  procedure  for  alias-free  simulations  of  compress¬ 
ible  turbulence.  The  main  difficulty  in  controlhng  aliasing  errors  stems  from  the 
division  by  density  to  yield  velocity  components  and  temperature  from  conservative 
variables.  If  the  equation  for  the  specific  volume  is  solved  instead  of  the  density 
equation  for  mass  conservation,  we  can  multiply  by  the  specific  volume  wherever 
division  by  density  is  needed.  The  equation  for  the  specific  volume,  v  =  1/p,  is: 


^  _  jdipuj) 

dt  ^  dxi 


(C.l) 


The  remaining  equations  are  the  momentum  and  energy  equations  in  conservative 
forms.  The  velocities  and  temperature  are  obtained  from  the  conservative  variables 
using 


Ui  =  V  (pu.) 


T  =  7v 


E  AP^k){p^k) 
^  2 


The  equation  of  state,  (3.17),  can  be  recast  as 


pv  = 


(t  - 

7 


(C.2) 


((7.3) 


Another  possible  source  of  alicising  errors  is  the  non-integer  power  law  depen¬ 
dence  of  viscosity  on  temperature.  The  effect  of  this  error  on  the  solution  remains 
to  be  investigated. 
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APPENDIX  D 

DIRECT  NUMERICAL  SIMULATIONS 
OF  SPATIALLY  EVOLVING  TURBULENCE 

Most  direct  or  large  eddy  simulations  of  turbulent  flows  have  been  performed 
with  periodic  boundary  conditions  [Rogallo  and  Moin  1984].  Direct  or  large  eddy 
simulations  of  flows  in  complex  geometries  require  “turbulent”  conditions  at  the 
inflow  boundary. 

This  Appendix  presents  a  method  of  generating  inflow  turbulence  with  a  pre¬ 
scribed  power  spectrum.  The  results  from  simulations  of  spatially  decaying  turbu¬ 
lence  are  then  compared  with  those  from  the  corresponding  temporal  simulation 
to  V2didate  the  method.  The  method  is  also  validated  by  comparison  with  the  ex¬ 
perimented  data  for  decaying  isotropic  turbulence  in  the  regime  where  the  energy 
spectrum  undergoes  a  self-similar  decay  [Ling  and  Huang  1970]. 

D.l  Method  of  Generating  Inflow  Turbulence 

For  simplicity,  we  consider  a  turbulent  flow  evolving  in  the  xj  direction,  the  mean 
flow  direction,  homogeneous  in  coordinates  X2,X3,  and  statistically  stationary  in 
time.  At  the  entrance  (ij  =  0)  of  the  computational  domain,  the  energy  spectrum 
of  a  flow  variable,  Eff,  is  prescribed  in  terms  of  frequency  and  two  transverse  wave 
numbers.  The  Fourier  coefficients,  are  prescribed  by  the  equation 

/(fc2,fc3,u>,t)  =  yjEff{k2,h,Lo)exp[i(f>r{k2,h,u;,t)],  (D.l) 

where  (f>r  is  the  phase  factor  and  i  =  y/-^.  The  functional  dependence  of  4>r  on  time, 
as  well  as  on  frequency  and  wave  numbers,  is  necessary  so  that  the  signal  generated 
is  not  periodic  in  time.  The  variable  t  in  f  denotes  an  element  of  the  ensemble  of 
realizations.  In  our  scheme,  (pr  is  changed  only  once  in  a  given  time  interval,  Tr,  at  a 
random  instance  by  a  random  amount,  A<f>r,  where  [  A<^r  I  is  bounded  by  a  prescribed 
value,  A(/>max-  The  randomized  temporal  dependence  of  each  phase,  <f>r{k2,k2,u),t), 
is  shown  schematically  in  Figure  D.l.  Because  the  phases  are  time-dependent,  the 
generated  signal  is  not  continuous  and  the  frequency  spectrum  of  the  turbulence 
signal  generated  only  approximates  the  target  spectrum,  Efj{k2,k2,u>).  The  fre¬ 
quency  and  amplitudes  of  the  random  phase  shifts  determine  the  smoothness  of 
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the  generated  signal.  If  this  temporal  dependence  is  weakened,  the  approximation 
of  the  target  spectrum  improves,  but  the  generated  signal  also  approaches  a  tem¬ 
porally  periodic  signal.  The  fluctuation  signal,  f{x2,X2yt),  is  obtadned  by  Fourier 
transforms  in  the  homogeneous  directions  (i.e.,  X2  and  x^)  followed  by  a  sum  over 
all  frequencies.  The  turbulence  signal  at  the  inflow  plajie  is  prescribed  by  adding 
the  fluctuation  signal,  fix2,X2,t),  to  the  mean  flow  profile,  F(x2,X3).  One  can 
easily  generalize  this  method  to  create  inhomogeneous  “turbulence”  as  well. 

The  approximation  of  the  target  spectrum  depends  on  the  choices  T*  —  u)oTr/2TT 
and  =  A4>max/2K,  where  Wo  is  the  frequency  of  the  peak  energy.  The  sen¬ 
sitivity  of  the  approximation  to  T*  and  shown  in  Figure  D.2.  The  target 

spectrum,  is  of  the  form  exp[— 2(a>/a;o)^].  As  the  temporal  dependence 

of  the  random  phase  increases  (small  T*  and  large  the  approximation  of 

the  target  spectrum  worsens.  The  region  of  disagreement  is  localized  in  the  non- 
energetic  frequencies.  In  the  high  cj  range,  the  spectrum  has  a  tail.  Figure 
D.3  shows  typical  turbulence  signals  generated  at  the  inflow. 

D.2  Simulation  of  Spatially  Decaying  Turbulence 

The  first  part  of  .s  section  demonstrates  the  performance  of  the  spatial  sim¬ 
ulation  and  compares  it  with  a  corresponding  temporal  simulation.  In  the  second 
part,  direct  numericcd  simulation  results  are  compared  with  the  experimental  data. 

D.2.1  Comparison  with  Temporal  Simulation  Results 

The  method  of  generating  inflow  turbulence  described  in  Section  II  is  used  to 
conduct  simulations  of  spatially  decaying  compressible  isotropic  turbulence.  The 
governing  equations  are  the  continuity  equation,  three  momentum  equations,  and 
the  energy  equation,  along  with  the  equation  of  state.  We  assume  the  fluid  to  be 
an  ideal  gas  (7  =  1.4)  with  zero  bulk  viscosity.  Viscosity  is  assumed  to  have  power- 
law  dependence  on  temperature,  ~  while  the  Prandtl  number,  Pr,  is 

kept  constant  at  0.70.  Special  attention  is  given  to  ensure  numerical  conservation  of 
mass,  kinetic  energy,  and  total  energy  in  the  inviscid  limit.  We  approximate  spatial 
derivatives  by  a  compact  finite-difference  scheme  [Lele  1990],  which  has  spectral- 
like  resolution  as  well  as  sixth-order  formal  accuracy.  Time  advancement  is  done 
explicitly  by  a  third-order  R\mge-Kutta  method.  Periodic  boundary  conditions 
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are  used  in  the  X2  and  13  directions.  The  mean  streamwise  velocity,  t/i,  is  kept 
supersonic  and  uniform  for  the  most  rigorous  application  of  the  outflow  boundary 
conditions.  In  the  example  presented,  the  mean  inflow  Mach  number,  M\  =  I/i/c, 
is  2.0  (c  is  the  speed  of  sound).  Simulations  conducted  for  subsonic  inflow  yield 
essentially  the  same  results.  Inflow  turbulence  is  generated  with  zero  density  and 
temperature  fluctuations,  with  a  three-dimensional  energy  spectrum  given  by 

f;(*:)~Jfc^exp[-2(ik/Jko)^].  {D.2) 


The  inflow  fluctuation  Mach  number  Mt  —  q/c  (q  =  where  is  a 

fluctuation  from  the  ensemble  averaged  velocity)  and  the  turbulence  Reynolds 
number  based  on  the  Taylor  microscale,  =  UoXi/u,  constitute  two  indepen¬ 
dent  parameters  of  the  simulation,  where  Uo  is  the  rms  fluctuation  in  a  velocity 
component.  We  consider  severed  cases  with  the  inflow  fluctuation  Mach  numbers 
Mi  =  0.519,  0.346,  0.173  and  Rex  =  25.0.  The  size  of  the  computational  domain 
is  27r  in  each  direction  with  64  grid  points.  Thus,  the  computational  wave  numbers 
are  integers,  and  we  use  the  energy  peak  wave  number,  ko  =  4,  in  prescribing  the 
inflow  spectrum.  For  comparison,  a  corresponding  computation  of  temporally  de¬ 
caying  turbulence  (Lee  et  al.  1991bj  is  conducted  with  the  initial  fluctuation  Mach 
number  and  Reynolds  number  Mt  =  0.346  and  Rex  =  25.0,  respectively.  Taylor’s 
hypothesis  is  used  to  convert  the  downstream  distance  from  the  inflow  boundary 
in  the  spatial  simulation  into  the  evolution  time,  i.e.,  t  =  x\ll]\. 

Figure  D.4  shows  evolution  of  velocity  derivative  skewness,  which  is  a  measure  of 
inertied  nordinearity  of  turbulence.  Skewness  varies  with  compressibility  as  well  as 
turbulence  Reynolds  number  [Tavoularis  et  al.  1978,  Erlebacher  et  al.  1990],  having 
a  value  of  about  —0.4  to  —0.6  for  isotropic  turbulence  at  Rex  ~  25.  Turbulence 
may  therefore  be  considered  realistic  beyond  time  t/rt  =  0.4,  where  the  turbulence 
time  scale  is  defined  as  =  Xjuo.  The  development  of  the  velocity  derivative 
skewness  for  the  spatially  decaying  turbulence  compares  favorably  with  that  for 
the  temporally  decaying  turbulence. 

Figures  D.5(a)  and  D.5(b)  present  comparisons  between  spatial  and  temporal 
simulations  for  one-dimensional  spectra  of  vorticity  and  dilatation  as  a  function 
of  the  transverse  wave  number  k2.  The  vorticity  spectra  agree  closely,  while  the 
dilatation  spectra  show  some  differences. 
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The  evolution  of  turbulent  kinetic  energy  is  shown  in  Figure  D.6.  The  turbu¬ 
lent  kinetic  energy  in  the  spatial  simulation  compares  favorably  with  that  in  the 
temporal  simulation.  The  same  kind  of  agreement  is  obtained  in  the  evolution  of 
rms  vorticity  and  turbulence  Reynolds  number.  These  agreements  contrast  with  a 
systematic  difference  in  the  statistics  of  dilatation.  Close  comparisons  are  generailly 
found  for  the  statistics  donrinated  by  the  incompressible  part  of  turbulence  [MoyaJ 
1952],  whereas  the  statistics  dominated  by  the  flow  compressibihty  tend  to  differ. 
The  level  of  rms  dilatation  of  the  spatial  simulation  when  Mt  =  0.346  is  lower  than 
the  corresponding  temporal  simulation  by  15  percent. 

The  deviation  in  the  dilatation  statistics  may  be  attributed  to  two  causes.  Firstly, 
disturbances  in  incompressible  turbulence  are  generally  advected  at  the  mean  flow 
speed,  Ui,  while  fluctuations  in  compressible  turbulence  are  converted  at  different 
speeds,  U\,Ui  -t-  c,  and  Vi  —  c.  Hence,  for  statistics  dominated  by  compressibil¬ 
ity  the  use  of  Taylor’s  hypothesis  may  be  inaccurate.  Secondly,  as  Figure  D.5(b) 
illustrates,  the  level  of  the  dilatation  spectrum  for  low  wave  numbers  is  higher  in 
the  temporally  decaying  flow  than  in  the  spatially  decaying  flow.  This  higher  level 
of  compressibihty  in  the  temporal  simulation  may  be  attributed  to  the  inability  of 
the  periodic  boundary  conditions  to  freely  radiate  the  acoustic  waves  generated  by 
turbulence.  The  existence  of  one  freely  radiating  boundary  (outflow  boundary)  in 
the  spatial  simulation  lowers  the  overall  intensity  of  acoustic  waves  trapped  in  the 
domain.  Because  of  this  difference,  caution  must  be  exercised  in  using  periodic  (or 
temporal)  simulation  databases  to  examine  compressibility-driven  quantities  such 
as  dilatation  dissipation  and  pressure  dilatation  correlation. 

D.2.2  Comparison  with  Experimental  Data 

Ling  and  Huang  [1970]  found  the  isotropic  turbulence  decay  between  microscale 
Reynolds  numbers  3  and  30  to  be  self-similar.  At  low  Reynolds  numbers  the  en¬ 
ergy  and  dissipation  scales  overlap  so  that  there  is  only  a  single  characteristic  leng'n 
scale.  Domaradzki  and  Rogallo  [1990]  confirmed  this  finding  through  direct  numer¬ 
ical  simulation  of  temporally  decaying  incompressible  turbulence  by  showing  that 
the  three-dimensional  energy  spectra,  E{k),  at  different  times  can  be  collapsed  us¬ 
ing  the  Taylor  microscale  as  the  characteristic  length  scale.  VVe  init’  Jized  a  direct 
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numericail  simulation  of  spatially  decaying  isotropic  turbulence  with  the  normalized 
inflow  energy  spectrum  given  by  Ling  and  Huang  [1970]  as 

=  aJb*(l  +  jb*)exp(-it*),  {D.3) 

where  k*  is  the  magnitude  of  the  normalized  wave  number,  k*  =  with 

k*  =  kiay/v{t  —  to),  a  —  3.162,  E*  —  E{k)y/u{i  —  to)/u'^^,  and  t  =  x\IU\  and  to 
are  the  decay  time  and  the  virtual  origin  of  the  decay  time,  respectively. 

The  corresponding  one-dimensional  energy  spectrum  is 

jEj  (fcj)  =  a  exp(  — fcj).  {DA) 

The  self-similarity  of  normahzed  spectra  in  {D.3)  and  {D.A)  is  based  on  the  as¬ 
sumption  that  turbulent  kinetic  energy  and  all  the  relevant  turbulence  length  scales 
evolve  like  and  respectively. 

In  the  simulation  of  the  experimental  conditions,  we  follow  the  same  numerical 
procedures  as  in  Section  D.2.1.  The  mean  streamwise  Mach  number,  fluctuation 
Mach  number,  and  Reynolds  number  at  the  inflow  are  Mi  =  1.20,M(  =  0.173,  and 
Re\  —  15.0,  respectively.  The  size  of  the  computational  domain  is  27r  in  each 
direction  with  64  grid  points,  and  the  energy  spectrum  peaks  at  k  =  A.  An  incoming 
fluid  particle  passes  through  the  computational  domain  in  time,  t  =  2.5t<,  where 
n  =  X/uo- 

Normalized  one-dimensional  spectra  at  different  downstream  positions  are  shown 
in  Figure  D.7.  Spectra  collapse  onto  the  experimental  results  with  some  deviation 
at  small  wave  numbers  because  of  the  limited  sample  size.  It  was  experimentally^ 
observed  that  turbulent  kinetic  energy  decays  like  (xj  —  Xo)~‘^,  i.e.. 


( 


X I  Xq 


{D.5) 


where  qo/2  is  the  turbulent  kinetic  energy  at  inflow  (i.e.,  x\  =  0)  and  Xq/Uiti  = 
—4.2  in  the  simulation.  Figure  D.8  compares  the  turbulent  kinetic  energy  decay 
in  the  simulation  with  that  expressed  by  {D.5).  The  evolution  of  turbulent  kinetic 


energy  in  the  numerical  simulation  is  in  excellent  agreement  with  the  experimental 
data. 

We  have  developed  a  method  of  generating  inflow  turbulence  fluctuations  for 
spatially  developing  turbulence  computations.  Using  this  method  we  performed 
direct  numerical  simulation  of  spatially  evolving  isotropic  turbulence.  The  com¬ 
puted  incompressible  turbulence  statistics  are  in  excellent  agreement  with  those 
from  the  corresponding  temporal  simulation  and  the  experimental  data.  The  al¬ 
gorithm  for  generating  inflow  turbulence  is  by  no  means  unique.  The  significant 
result  is  that  one  apparently  can  compute  spatially  developing  turbulence  with  the 
accuracy  typical  of  present  temporally  evolving  simulations. 
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different  wave  numbers  and  frequencies. 


Figure  D.2.  Approximation  of  the  spectrum  by  signals  with  time-dependent 
random  phases:  - target  spectrum, - T*  =  =  1/20, - T*  = 


2)  ^<^mcLX  —  l/lOi  T*  —  —  1/20. 
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Figure  D.7.  Comparison  of  the  normalized  energy  spectrum  of  Ling  and  Huang 
[1970]  given  by  {DA)  (solid  line)  and  the  energy  spectra  obtained  in  the  simulation 
at  five  different  positions  (symbols):  o  ,  xj  =  0;  A  ,  xj  =  7r/2;  +  ,  xj  =  tt;  x  , 
xj  =  37r/2,  and  o  ,  xj  =  27r. 
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APPENDIX  E 


PARAMETER  LIMITATIONS  FOR  DIRECT  NUMERICAL  SIMULATION 
OF  SHOCK  TURBULENCE  INTERACTION 


In  our  direct  numerical  simulations  of  shock  turbulence  interaction,  we  have  re¬ 
quired  the  resolution  of  all  the  relevant  scales  of  turbulence  and  the  shock  wave.  The 
resolution  requirement  for  incompressible  turbulence  is  far  better  established  than 
that  of  compressible  turbulence.  For  weakly  compressible  turbulence  with  Mt  <  0.2 
where  there  are  no  significant  compressibility  effects  on  turbulence,  resolution  re¬ 
quirements  for  compressible  turbulence  is  comparable  to  that  of  incompressible 
turbulence.  Some  simulations  of  compressible  turbulence  at  very  high  Mt  [Kida  et 
al.  1990,  Blaisdell  et  al.  1990,  Sarkar  et  al.  1991]  were  not  successful  in  resolving 
changes  of  the  flow  variables  across  the  eddy  shocklets.  Since  their  main  interests 
were  to  study  the  effect  of  compressibility  on  turbulence  evolution,  resolving  the 
eddy  shocklets  may  not  have  been  a  critical  factor.  In  this  study  we  investigated 
the  effect  of  turbulence  on  the  shock  wave  structure  as  well  as  the  effect  of  a  shock 
wave  on  turbulence.  Therefore,  proper  resolution  of  the  shock  wave  is  necessary 
for  accurate  results.  Limitations  on  the  availability  of  computer  resources  leads  to 
limitations  on  the  choice  of  physical  parameters  in  the  simulation,  and 

Mt. 


E.l  Requirements  for  Resolution  and  Sample  Size 

As  shown  in  Section  3.3.1,  proper  resolution  of  the  shock  wave  structure  requires 
a  grid  distribution  which  places  at  least  three  points  inside  the  shock  wave  whose 
thickness  6s  is  defined  by 

6 

Thus,  6s  >  3(Ai|)n^in,  where  (Aii)n,in  is  the  typical  streamwise  grid  spacing  near 
the  shock  wave.  For  weak  shock  waves,  the  shock  thickness  6s  is  estimated  as 


8  v' 

3c{M  -  1) 


(E.l) 
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[Thompson  1984],  where  u  and  c  are  the  upstream  kinematic  viscosity  and  speed 
of  sound,  and  M  is  the  maximum  instajitaneous  upstream  Mach  number,  M  cr 
4-  Mt,  where  Mt  is  the  fluctuation  Mach  number  at  the  inflow. 

In  the  numerical  simulation,  artificial  turbulence  generated  at  the  inflow  evolves 
into  “realistic”  turbulence  as  it  approaches  the  shock  wave.  In  Appendix  D,  evolu¬ 
tion  of  artificial  inflow  turbulence  into  realistic  turbulence  was  found  to  take  place 
in  a  distance  which  determines  the  minimum  streamwise  computational  box  size 
upstream  of  the  shock  wave.  Using  the  results  of  Appendix  D  and  applying  Taylor’s 
hypothesis,  we  get  the  following  relation: 


t 

rt 


A/uf) 


>  0.4, 


(£.2) 


where  A  is  the  longitudinal  Taylor  microscale,  and  rj  =  X/uq  is  a  turbulence  time 
scale.  Therefore,  the  minimum  streamwise  computational  box  size  upstream  of  the 
shock  wave  is 


U_  0-4A 
^  Uo/Ul' 


(E.Z) 


where  Uo/Ui  is  the  streamwise  turbulence  intensity  at  the  inflow.  The  streamwise 
distance  downstream  of  the  shock  wave,  was  chosen  to  be  the  same  as  .  The 
streamwise  computational  box  size  Tj  can  now  be  expressed  as 


Li  = 


£)  0.8A 

^  "  uo/Ui  ■ 


This  relation  can  be  rewritten  as 


_  0.46A 
‘  ^  Mt/Mf' 

using  the  definitions  of  Mt  and  M['^ 


(EA) 


Mt  = 
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In  order  to  eliminate  the  contamination  of  long-time-averaged  statistics  across 
the  shock,  the  mean  shock  front  should  not  drift  in  space.  As  the  upstream  fluc¬ 
tuation  Mach  number,  Mt,  becomes  comparable  to  —  1,  the  mean  shock  drift 
speed  increases.  Therefore,  the  fluctuation  Mach  number  should  be  bounded  by 
Mf  -  1. 

Compressible  turbulence  of  Mt  <  0.2  does  not  contain  eddy  shocklets  [Lee  at  al. 
1991b]  and  can  be  resolved  by  a  grid  similar  to  that  in  incompressible  turbulence  at 
the  same  Reynolds  number.  In  our  simulations  turbulence  Mach  number  is  always 
lower  than  0.2  except  possibly  just  behind  the  shock  wave. 

The  size  of  the  computational  domain  and  the  number  of  grid  points  in  transverse 
directions  are  determined  to  satisfy  the  requirements  of  sample  size  and  resolution 
of  turbulence.  Based  on  the  form  of  inflow  turbulence  spectrum  given  in  (3.32), 
the  integral  scale  1^22, 2  defined  in  (4.7),  the  longitudinal  Taylor  microscale  A,  and 
the  Kolmogorov  length  scale  are,  respectively, 

^0^22,2  — 

koX  =  2  {E.5) 


where  Rex  is  Taylor  microscale  Reynolds  number  defined  in  (4.1). 

To  have  a  sufficient  sample  of  large  scale  structures,  the  computational  box  in 
the  transverse  directions  were  chosen  to  be  larger  than  about  ten  integral  scales, 

^2  —  ^3  ^  10A22,2-  iE.6) 

To  resolve  the  smallest  scale  of  turbulence  which  is  comparable  to  rjj^,  we  chose 
the  mesh  size  in  the  transverse  direction  as 

Ai2  =  Ax3  <  2Tjf^ 

[Lee  an  Reynolds  1985],  or  using  {E.5) 
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jfcoAx2  =  0.508/2e;^ {E.l) 

Using  the  expression  for  in  {E.o),  the  ratio  of  the  cutoff  wave  number  kc  to  the 
energy  peak  wave  number  kg  can  be  expressed  as 

where  the  relation 

kc  =  ^N2  {E.9) 

L2 

is  used  {Nn  is  the  number  of  grid  points  in  the  X2  direction).  Combining  this  with 
(£’.5),  (£.6),  (£.8)  and  (£.9),  the  number  of  grid  points  in  the  transverse  directions 
can  be  expressed  as 


N2  =  N3>  12£ey^  (£.10) 

The  number  of  grid  points  in  the  transverse  directions  according  to  (£.10)  satisfies 
the  requirements  of  sample  size  and  fine  scale  resolution. 

A  non-uniform  mesh  is  used  in  the  streamwise  direction.  The  minimum  mesh 
size  to  resolve  the  shock  wave  is  obtained  from  (£-1)  as 


^o(Aii)n,in  — 


16  Mj  1 

9^3  Rex  Mf  +  Mt-l' 


The  maximum  streamwise  mesh  spacing  is  equal  to  the  mesh  spacing  in  the  trans¬ 
verse  directions. 

Turbulence  statistics  are  obtained  with  averaging  over  a  time  i  such  that 


Uot 

X 


>4, 


and  the  maximum  time  step  of  the  simulation  is  determined  by  the  numerical 
stability  condition  of 


CFL  =  XX’ 

(  )inin 
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where  CFL  is  the  Courant-Friedrichs-Lewy  number.  Since  the  smallest  grid  spacing 
is  in  the  ij-direction,  the  total  number  of  time  steps  required  in  each  simulation 
Nf  is  estimated  to  be 


t  ^  4V3  A  +  1 

At  -  CFL  (AxiU„  Mt 


(EM) 


By  using  6s  >  3(Axi)jj,in,(^-l),(£-5),(£.8)  and  {E.9),  one  can  rewrite  (£^.11)  as 

W,  ^(Mf +  M,-1)(M['  +  1).  (£.12) 


In  the  present  study  with  the  choice  of  parameters  given  in  Table  4.1,  the  number 
of  time  steps  ranges  from  10,000  to  30,000. 

The  separation  between  the  typical  turbulence  length  scale  and  the  shock  wave 
thickness  requires  a  large  grid  stretching  in  the  streamwise  direction.  Furthermore, 
acoustic  waves  should  be  accurately  resolved  in  the  region  occupied  by  the  shock 
wave,  which  excludes  the  possibility  of  using  an  implicit  time  advancement  to  take 
a  large  time  increment.  Explicit  time  advancement  with  a  small  time  increment  is 
the  main  reason  for  the  large  CPU  time  required  for  direct  numerical  simulation  of 
reahstic  shock/turbulence  interaction. 

In  a  three-dimensional  direct  numerical  simulation,  storage  of  all  the  flow  vari¬ 
ables  in  the  core  requires  a  large  memory.  On  the  Cray  Y-MP/832,  this  problem  is 
alleviated  by  using  the  Solid  State  Device  (SSD)  I/O.  In  the  present  computations, 
the  SSD  I/O  takes  about  1.5%  of  the  total  CPU  time  used. 

The  code  performance  is  about  130  megaflops  with  129  x  64  x  64  grid  points,  and 
uses  about  45  CPU  seconds  of  Cray  Y-MP/832  for  each  time  step.  It  takes  about 
100  to  200  CPU  hours  of  Cray  Y-MP/832  to  obtain  adequate  turbulence  statistics 
for  one  case. 


E.2  Limitations  of  the  Physical  Parameters:  Rex,  MW  Aft 

To  dii^^ass  the  limitations  of  the  physical  parameters,  we  use  as  a  reference  the 
simulated  case:  Re^  =  25.0,  =  1.20,  and  Aft  —  0.173.  In  realistic  problems 
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with  shock/turbulence  interaction,  the  turbulence  Reynolds  number  is  higher,  and 
the  shock  strength  is  usually  stronger  with  a  wider  range  of  Mt- 

Increase  in  the  Reynolds  number,  Rej^,  results  in  an  increase  in  the  number  of 
grid  points  (ref.  (£.10))  in  all  directions  which  leads  to  an  increase  in  the  CPU  time 
per  time  step.  Increasing  the  upstream  Mach  number,  Alf ,  requires  a  larger  grid 
stretching  ratio  in  the  streamwise  direction  due  to  the  increased  scale  separation 
between  turbulence  and  the  shock  wave  (ref.  (£.1))  and  more  CPU  time  per  time 
step.  Increasing  the  Mach  number  to  M[^  =  1.40  requires  at  least  twice  the  CPU 
time  required  for  the  reference  case. 

Lowering  the  upstream  fluctuation  Mach  number,  A/^,  with  fixed  upstream  tur¬ 
bulence  length  scale  requires  a  longer  computational  box  size  in  the  streamwise 
direction  (ref.  (£.4)),  which  leads  to  the  increase  of  the  number  of  grid  points. 
Higher  Mt  requires  more  grid  points  in  the  streamwise  direction  to  resolve  the  in¬ 
stantaneous  shock  wave  structure  (ref.  (£.1)).  The  upper  limit  of  Mt  is  usually  set 
by  the  condition  of  shock  stationarity,  Mt  <  M^  -  1. 

As  shown  in  Table  4.1,  the  ranges  of  inflow  parameters  used  in  the  simulations 
were 


16  <Rex  <  25, 
1.05  <A/f  <  1.20, 
0.087  <Mt  <  0.173. 
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APPENDIX  F 


TURBULENCE  STATISTICS  IN  KINEMATIC  OSCILLATION 
OF  A  PLANE  SHOCK  WAVE 

In  the  interaction  of  turbulence  with  a  shock  wave,  the  shock  front  is  distorted 
and  undergoes  an  oscillatory  movement.  The  oscillation  of  the  shock  front  produces 
an  intermittent  time  history  of  the  flow  variables  at  a  fixed  point  in  a  reference 
frame  fixed  at  the  mean  shock  position,  and  it  leads  to  overprediction  of  turbulence 
statistics  (see  Section  4.1  and  [Debieve  et  al.  1986]).  This  large  fluctuation  (which  is 
not  turbulence)  is  driven  by  the  upstream  turbulence,  and  undergoes  rapid  viscous 
decay  since  its  time  and  length  scales  are  small. 

In  order  to  identify  the  effects  of  the  shock  oscillation  on  turbulence  statistics,  a 
plane  shock  wave  was  moved  back  and  pro  in  a  sinusoidal  fashion  in  the  streamwise 
direction,  and  the  budget  of  terms  in  the  equation  in  (4.9)  was  computed.  The 
profiles  of  the  flow  variables  across  the  shock  wave  were  obtained  as  the  laminar 
solution  of  the  one-dimensional  Navier-Stokes  equations  for  =  1.20  (see  Section 
3.2.2  for  more  details  on  the  solution  procedure).  And  the  solution  U  =  (/>,ui,p)^  is 
expressed  in  terms  of  the  relative  position  with  respect  to  the  shock  center  position 
Xs  as  (ij  —  15)/^,,  where  83  is  the  shock  thickness.  The  oscillatory  movement  of 
the  shock  wave  is  emulated  by  externally  driving  the  shock  wave  to  move  back  and 
forth  in  time  as 


is(t)  =  a:s(0) -i- Qj(5,/(t),  (F.l) 

where  0363  is  the  spatial  amplitude  of  the  oscillation,  and  f{t)  is  a  periodic  function 
with  f{t)  ~  0.  Statistical  samples  were  taken  at  fixed  points  which  were  apart 
by  about  a  tenth  of  the  shock  thickness.  The  flow  variables  at  a  sampling  point  near 
the  shock  wave  vary  in  time  due  to  changes  in  the  relative  positions  with  respect  to 
the  shock  wave.  The  sampling  points  far  away  from  the  shock  wave  are  not  affected 
by  the  shock  wave  oscillation.  The  sampling  time  interval  was  chosen  so  that  100 
samples  were  taken  per  oscillation  period  of  the  shock  wave.  The  statistics  were 
not  found  to  be  sensitive  to  the  choice  of  the  parameter,  Qj,  and  the  function,  f{t). 
In  the  following,  a  sinusoidal  function  is  used  for  /(/)  with  —  0.25. 


197 


Figure  F.l  shows  the  statistics  of  the  streamwise  velocity,  dilatation,  pressure, 
density,  and  temperature.  As  found  in  Section  4.2,  all  the  statistics  peak  inside  the 
zone  of  shock  oscillation.  The  peak  values  were  also  comparable  to  those  from  the 
direct  numerical  simulation  (see  Figures  4.4  and  4.24). 

Figure  F.2  shows  the  budget  of  terms  in  the  equation.  The  behaviors  of  all 
the  terms  in  the  equation  were  found  to  be  consistent  with  those  computed  from  the 
direct  numerical  simulations  (see  Figure  4.8(a)).  The  pressure  work  term, 
was  found  to  be  the  dominant  term  inside  the  shock  wave. 

Figure  F.3  shows  the  decomposition  of  the  pressure  work  term  into  the  pressure 
transport  term,  and  the  pressure-dilatation  correlation,  p'u'^  Again, 

behaviors  of  the  two  decomposed  terms  are  in  agreement  with  the  results  from  the 
simulations  (s<“e  Figure  4.7(a)). 

It  can  be  concluded  that  the  large  levels  of  fluctuations  in  the  shock  zone  in 
the  direct  numerical  simulations  were  mcunly  due  to  the  oscillatory  movement  of  a 
plajie  shock  wave. 
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Figure  F.l.  Effects  of  kinematic  shock  wave  oscillation  on  turbulence  statistics 
Prms,  3rms- 
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Figure  F.3.  Effects  of  kinematic  shock  wave  oscillation  on  the  pressure-work  term  decomposition: 
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APPENDIX  G 


EFFECT  OF  A  REFINED  OUTFLOW  BOUNDARY  CONDITION 


As  described  in  Section  3.2.4,  Thompson’s  non-reflecting  boundary  condition 
[1987]  was  used  in  the  present  work.  This  boundary  condition  was  successful  in 
suppressing  the  reflection  of  nonphysical  acoustic  waves  in  test  problems,  where 
the  vortical  and  entropy  waves  were  passed  through  the  outflow  boundary  (see  Sec¬ 
tion  3.3.2).  However,  in  the  simulations  of  shock  wave  turbulence  interaction,  the 
statistics  which  are  associated  with  acoustic  waves,  such  as  pressure  work  (Figure 
4.6(b))  and  dilatation  (Figure  4.10),  show  anomalous  behaviors  near  the  outflow 
boundary.  In  order  to  assess  the  extent  of  influence  of  downstream  boundary  con¬ 
ditions  on  turbulence  evolution,  the  more  refined  boundary  conditions  of  Giles 
[1990]  were  implemented  and  turbulence  statistics  from  these  computations  were 
compared  with  those  using  Thompson’s  boundary  conditions.  In  this  Appendix,  a 
brief  description  of  Giles’  boundary  condition  is  given  first,  followed  by  a  compari¬ 
son  of  the  statistics  using  Giles’  boundary  condition  with  those  using  Thompson’s 
boundary  condition. 

G.l  Description  of  Giles’  Boundary  Condition 


The  original  derivation  of  Giles’  boundary  condition  is  based  on  the  analysis  of 
the  linearized  Euler  equations;  here  the  viscous  terms  were  added.  We  begin  with 
Giles’  analysis  on  the  three-dimensional  Euler  equations  which  can  be  written  in 
terms  of  the  primitive  variables  as 


dU  ,dU  „  dU 
dt  ^dxi  ^dx2 


+  c 


dU 

dxz 


-0, 


where  U  =  {6p,8ui,6u2,6u:i,6p)'^ , 
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The  elements  of  the  vector  U  represent  perturbations  from  uniform  flow  condi¬ 
tions,  and  the  matrices  A,  B,  and  C  are  evaJuated  using  uniform  flow  conditions. 
We  consider  wave-like  solutions  of  the  form 


U{xi,X2,X2,t)  =  exp[z(fcia:i  -f-  k2X2  -f-  ^313  -  u>t)]u^,  {G.5) 

where  is  a  constant  column  vector.  Substituting  this  into  the  differential  equa¬ 
tion  (G.l),  we  obtain 


{-u^I +  kiA  +  k2B +  k2C)u^  =  0,  (G.6) 

which  has  a  nontrivicd  solution,  provided  that 

det(— u;/  -f-  kiA  +  k2B  -I-  k^C)  =  0.  [G-^) 

The  vector  is  cilso  an  eigenvector  of  the  matrix 

H  =  A-\~u;I  +  kiA  +  k2B  +  k:iC)  (G.8) 

corresponding  to  the  eigenvalue  fcj. 

Suppose  that  the  differential  equation  is  to  be  solved  in  the  domain  xj  <  Li, 
and  one  wants  to  construct  boundary  conditions  at  xj  =  to  minimize  or  ideally 
prevent  the  reflection  of  outgoing  waves.  At  the  boundary  X]  =  Lj,  f7  can  be 
decomposed  into  a  sum  of  Fourier  modes  with  different  values  of  ^25 ^3) 
Consider  a  single  wave  with  particular  choice  of  fc2>^3>  ^^^d  w.  In  this  case,  the 
most  general  form  for  U  is 
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[7(Li,Z2,*3>0  = 


t(fcj*2  +  fc3*3-w0 


(G.9) 


L  n 

where  (ki)n  is  one  of  the  roots  of  the  dispersion  relation  for  the  given  values  of 
k2,k^,  and  w,  and  is  the  corresponding  right  eigenvector. 

The  ideal  nonreflecting  boundary  condition  would  specify  Cn  =  0  for  each  n  that 
corresponds  to  an  incoming  wave.  The  construction  of  such  a  boundary  condition 
requires  the  vector  which  is  an  eigenvector  of  .  It  is  well  known  from  linear 
algebra  that  the  eigenvectors  of  H  and  corresponding  to  different  eigenvalues 
are  orthogonal,  that  is, 


=  0,  (G.IO) 

where  and  are  the  left  and  right  eigenvectors  corresponding  to  different 
solutions  {ki)m  and  (ki)n  of  the  dispersion  relation,  (G.7). 

At  the  boundary,  ij  =  LI,  orthogonality  leads  to 


(vifu  =  w) 


L  rn 


RJ(^i)mLi 


On 


“?]  e” 


(kl)nLi  i(k2X2  +  k3X3-U/t) 


(G.ll) 


Therefore,  an  equivalent  specification  of  nonreflecting  boundary  condition  is 


(vifu  =  0  (G.12) 

for  each  n  corresponding  to  an  incoming  mode. 

In  principle,  these  exact  boundary  conditions  can  be  implemented  in  a  numerical 
method.  The  problem  is  that,  in  general,  depends  on  I2  =  ^2/^^  h  =  k^/u, 
and  the  implementation  would  involve  Fourier  transforms  in  Z2  *3  Laplace 
transform  in  t.  Computationally,  this  is  both  difficult  and  expensive  to  implement. 
In  the  following,  an  approximation  used  for  general  situations  is  described. 
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A  sequence  of  approximate,  nonreflecting  boundary  conditions  can  be  obtained 
by  expanding  in  a  Taylor  series  a^  a  function  of  I2  and  as 


V^ih.h)  =  0^(0, 0)  +  h^iO,Q)  +  <3^(0, 0)  +  0(llhh,ll).  (G.13) 


The  first-order  approximation,  obtained  by  keeping  only  the  leading  term,  gives 
Thompson’s  boundary  condition.  The  second-order  approximation  is. 


*^^<(0,0)  +  hM(0,0)' 


U>  5/2 


o)  9/3 


U  =  0. 


(G.14) 


Multiplying  by  —iu)  and  replacing  1^23*^33  9/8x2, d/dx^,  and  —  d/dt, 

respectively,  gives 


(»if(o.o)) 


dt 


dh 


8x2 


((7.15) 


This  is  a  local  boundary  condition  of  the  same  differential  order  as  the  governing 
equations.  These  boundary  conditions  are  only  approximately  nonreflecting  and 
may  produce  nonphysical  reflections  of  outgoing  waves  for  which  I2  and/or  are 
far  from  zero. 

Using  ((7.2)-(G'.4),  the  dispersion  relation  ((7.7)  for  the  system  of  equations, 
((7.1),  can  be  written  as 


(uiki  -f-  U2^2  +  “3^3  —  |^(wi^l  +  +  “3^3  ~  '^)^  ~  +  ^2  'f"  ^1) 

Three  of  the  five  roots  are  identical,  i.e.. 


=  0. 
(G.16) 


(*^1)1, 2, 3  = 


u>  —  162^2  ~  ^3^3 
uj 


((7.17) 
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(G.18) 

(G.19) 


where 

D  -  -  (c2  -  u|)(fc|  +  kl)/{u}  -  U2k2  -  {G.20) 

For  0  <  uj  <  c,  which  corresponds  to  subsonic  flow  normal  to  the  boundary,  the 
fifth  root  is  a  left-travelling  wave,  provided  the  correct  branch  of  the  complex  square 
root  function  is  used  in  defining  D.  {D  is  defined  as  the  positive  root  for  >  0 
and  as  the  one  with  negative  imaginary  part  otherwise,  so  that  it  represents  a 
left-traveling  wave  with  finite  amplitude.)  This  wave  is  of  interest  in  implementing 
Giles’  boundary  condition  at  the  subsonic  outflow  boundary.  The  left  null-vector 
defined  in  (G.9)  for  (^1)5,  which  correspond  to  an  upstream  travelling  pressure 
wave,  is 


T  'i 

(^5  )  =  {0,pc{-u}  +  U2k2  +  U3k3),-pcuik2,-pcuik3,{u;-U2k2-U3k:i)J'>).  (G.21) 
The  second-order  approximation  of  (G.21)  in  the  form  of  (G.15)  is 


d\J 

iO,~pc, 0,0,1)—  -(0,pcu2,-pcui,0,-U2)^  -  (0,pcu3,0, -pcui, -U3)—  =  0. 

(G.22) 

For  convenience  of  implementation  and  for  comparison  with  Thompson’s  bound¬ 
ary  condition,  we  define  one-dimensional  characteristic  variables, 


C  =  (ci,C2,C3,C4,C5)^  , 
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and 


/ 

-c2 

0 

0 

0 

A 

0 

0 

pc 

0 

0 

u 

c  = 

0 

0 

0 

pc 

0 

0 

pc 

0 

0 

1 
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0 

-pc 

0 

0 

ly 

(  -l/c2 

0 

0 

1/  2c2 

1/  2c2)  a 

0 

0 

0 

l/(2pc) 

-l/(2pc) 

0 

l/(pc) 

0 

0 

0 

c 

0 

0 

i/M 

0 

0 

\  0 

0 

0 

1/2 

1/2  j 

(G.23) 


(G.24) 


We  can  express  the  boundary  condition  ((7.22)  using  a  characteristic  variable  as 


dc^ 

dt 


+  (0,ui,0,0,U2) 


dc 

dx2 


+  (0,0,ui,0,U3) 


dc 

8x3 


=  0. 


((7.25) 


The  corresponding  expression  for  Thompson’s  boundary  condition  is 


^  O 

+  (0,c,0,0,U2)~  +(0,0,c,0,U3)^ 


=  0, 


((7.26) 


which  simply  replaces  the  velocity  uj  in  ((7.25)  with  c. 

The  actual  boundary  conditions  used  in  the  present  code  were  those  given  by 
(6'.22)  with  viscous  terms  added.  The  implementation  of  Thompson’s  boundary 
conditions  also  included  the  viscous  terms. 


G.2  Comparison  of  Turbulence  Statistics 


The  code  used  in  the  present  study  was  modified  to  include  the  more  accurate 
Giles’  boundary  condition  instead  of  Thompson’s  boundary  condition.  One  of  the 
saved  fields  for  the  case  A  which  used  Thompson’s  boundary  condition  was  used 
as  the  initial  field.  Except  for  the  outflow  boundary  condition,  the  same  algorithm 
was  used.  After  a  transient  period,  statistic2d  samples  of  flow  variables  were  ac¬ 
cumulated.  Turbulence  statistics  are  compared  with  those  from  case  A  which  are 
computed  from  samples  taken  for  exactly  the  same  time  interval. 
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Figure  G.l  compares  the  evolution  of  mean  square  vorticity  for  the  two  different 
outflow  boundary  conditions.  The  difference  is  negligible  throughout  the  domain 
which  indicates  that  the  improved  boundary  conditions  for  removing  acoustic  wave 
reflections  have  virtually  no  effect  on  the  evolution  of  vortical  waves. 

Figure  G.2  compares  the  evolution  of  the  streamwise  velocity  fluctuations.  The 
anomcdous  increase  towards  the  outflow  boundary  is  reduced  by  applying  Giles’ 
boundary  condition.  This  improved  behavior  is,  however,  localized  near  the  outflow 
boundary  and  does  not  affect  the  evolution  of  turbulence  downstream  of  the  shock. 

Apparently,  the  anomalous  behavior  of  the  streamwise  velocity  fluctuations  near 
the  outflow  boundary  is  caused  by  the  pressure  work  term  (see  Figure  4.8(a)).  Fig¬ 
ure  G.3  compares  the  statistics  of  pressure  work  term,  —p'-u'!  in  the  TKE  trans- 
port  equation.  The  sudden  increase  in  the  pressure  work  term  towards  the  outflow 
boundary  is  reduced  by  half  by  using  Giles’  boundary  condition  instead  of  Thomp¬ 
son’s  boundary  condition.  This  improvement  is  also  localized  near  the  outflow 
boundary. 

As  shown  in  Figure  4.10,  the  variance  of  the  fluctuating  dilatation  increases  by 
a  factor  of  10  near  the  outflow  boundary  due  to  nonphysical  reflection  of  acoustic 
waves.  Figure  G.4  shows  the  improvement  resulting  from  removing  the  acoustic 
wave  reflection  by  the  use  of  Giles’  boundary  condition. 

In  conclusion,  Thompson’s  boundary  condition  used  in  the  present  study  is  found 
to  be  acceptable  because  reflections  of  acouttic  waves  at  the  outflow  boundary  do 
not  appear  to  affect  turbulence  statistics  downstream  of  the  shock. 
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APPENDIX  H 


DRIFT  IN  THE  SHOCK  POSITION  AND  THE  OUTFLOW  CONDITION 


In  our  numerical  setup,  there  is  no  external  mechanism  to  fix  the  position  of 
the  shock  wave.  From  numerical  experiments  with  a  one  dimensional  shock  wave 
interacting  with  a  sinusoidal  en'ropy  wave,  we  found  that  poor  resolution  of  the 
shock  and  a  high  amplitude  fluctuation  lead  to  the  shock  movement.  If  the  shock 
wave  drifts  in  the  streamwise  direction,  the  inflow  mean  Mach  number  is  not  the 
true  upstream  Mach  number  of  the  system.  Moreover,  the  statistics  obtained  at 
a  fixed  point  in  space  become  contaminated  because  the  relative  distance  to  the 
shock  wave  is  changing  in  time.  The  drift  of  the  shock  wave  position  in  time  is 
represented  as  the  mismatch  of  the  upstream  and  downstream  mass  fluxes  in  the 
chosen  reference  frame  of  the  simulation,  which  reflects  the  fact  that  the  mean 
turbulent  shock  propagation  speed  is  different  from  the  specified  laminar  shock 
propagation  speed.  The  shock  drift  speed  can  be  related  to  the  the  mass  flow  rate 
difference  by  integrating  the  continuity  equation  in  time  and  space,  as 


III  [p(<)  -  p(<o)l  dx\dx2dx-i  -  -  an  {pui)^  —  dx2dx^dt. 

{HA) 

The  left  side  of  {HA)  can  be  expressed  in  terms  of  the  mean  shock  drift  speed  Us 

1 1  J  ~  nu)]dxidx2dx:i  -{p^ -p^n^'sL2L2{t  -  to),  {H.2) 


as 


which  represents  the  mass  decrease  in  the  computational  domain  due  to  the  mean 
drift  of  the  shock  wave  (The  shock  wave  drift  in  the  +X]  direction  leads  to  a  decrease 
of  the  total  mass  in  the  computational  domain).  The  statistically  averaged  quantity 
~p  is  defined  as 

It[  1 1  P^^2dx3dt 
^  -  to) 
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The  right  side  of  (^.1)  can  be  rewritten  in  terms  of  the  statisticcdly  averaged 
quantities  as 


{pui)^  dx2dx^dt  =  —{pui 


D 


^^)L2L3{t-to).  {H.3) 


Using  {H.))  -  (i/.3),  the  mean  shock  drift  speed  can  be  obtained  as 


Us 


pui 


D 


pui 


U 


-pD^-pU 


{HA) 


Figure  H.l  shows  the  mismatch  in  the  mass  fluxes  across  the  shock  wave  for  the 
worst  case  among  the  simulations.  The  shock  drift  speed  in  this  case  is  about  0.7% 
of  the  average  upstream  speed. 

Since  the  flow  variables  at  the  outflow  boundary  are  not  explicitly  specified,  they 
may  drift  in  time.  For  the  simulations  of  subsonic  inflow-outflow  condition,  Poinsot 
et  al.  [1990]  proposed  to  add  an  extra  term  to  Thompson’s  boundary  condition  in 
order  to  force  the  exit  pressure  to  relax  to  an  equilibrium  pressure.  Figure  H.2 
shows  the  time  history  of  the  mean  exit  pressure  over  three  eddy  turnover  times. 
The  drift  in  the  mean  exit  pressure  is  about  0.2%  of  the  mean  exit  pressure,  which 
is  not  significant.  Therefore,  no  modifications  were  introduced  to  the  Thompson 
boundary  condition  at  the  outflow. 
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Figure  H.l.  Mass  flowrate  across  the  shock  wave  for  case  A 
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